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SUMMARY 


This  report  presents  the  current  status  of  work  under  ARPA  order  2089, 
in  the  weather  models,  with  a  view  toward  climatological  application. 

Vector  Spherical  Harmonics  are  reviewed  in  the  context  of  solving 
partial  differential  equations  in  epherical  coordinates.  The  has  c 
equations  governing  global  weather  models  are  presented.  Vector  Spherical. 
Harmonics  are  applied  to  these  equations  in  i spherical. ‘ ^  lar 

is  shown  how  separation  of  variables  is  achieved  in  the  sense  that  angular 

dependence  is  eliminated. 

Angular  momentum  is  expressed  in  terms  of  Vector  Spherical  Harmonics. 
The  equation  for  the  conservation  of  angular  momentum  takes  on  a  very 
simple  and  elegant  form  in  this  representation. 

.  The  changes  of  variable  necessary  to  introduce  pressure 
and  tf-coordinates  are  carried  out  and  the  consequences  of  these  changes 
discussed.  The  results  are  used  to  formulate  the  Mintz-Arakawa  model  in 
terms  of  Vector  Spherical  Harmonics  in  detail. 

The  transform  method  which  will  be  used  in  the  numerical  solution  rf 
the  model  is  reviewed. 

The  functional  specifications  for  the  computer  programs  to  implement 
the  above  results  are  included,  as  well  as  a  discussion  on  the  proposed 
implementation  methodology • 

How  this  transform  method  applies  to  Vector  Spherical  Harmonics  is 
shown. 


Introduction 


This  report  presents  the  main  results  of  the  first  six  months  work  supported 
by  ARPA  Contract  No.  DAHC  04-72- C -0020  on  global  weather  modeling  using 
vector  spherical  harmonics  (vsh).  The  work  to  date  has  been  mainly  mathe¬ 
matical.  The  primitive  equations  have  been  formulated  in  terms  of  vsh. 

The  non-linear  terms  can  all  be  expressed  in  terms  of  sums  of  products  of 
various  vector  coupling  coefficients.  These  are  being  developed,  but  are 
not  reported  on  here.  The  transform  method  to  be  used  for  the  computation 
does  not  require  the  computation  of  product  terms  in  vsh.  Thus  the  formu¬ 
lation  in  this  report  suffices  for  the  transform  approach.  However,  these 
explicit  expansions  are  important  for  checking  the  transform  method  results 
and  for  analytic  work.  We  expect  to  have  this  formulation  ready  for  the  nev+ 
six  month  report. 

The  Mintz-Arakawa  model  has  been  formulated  in  terms  of  <£  coordinates 
and  a  tentative  approach  using  powers  of  has  been  developed.  Work  on 
correlating  actual  data  with  vsh  is  in  progress.  The  nature  of  the  radial 
dependence  of  the  spectral  coefficients  is  being  investigated  in  this  way. 

A  formulation  using  vsh  with  pressure  coordinates  is  presented.  A  brief 
discussion  of  work  on  angular  momentum  in  progress  is  also  shown. 

Although  the  program  design  for  the  computer  calculations  has  been  done 
with  other  funds  and  is  not  necessary  for  fulfilling  tlie  contract  commitments, 
we  thought  it  of  sufficient  interest  to  be  included  in  this  report.  It  will  enhance 
the  work  on  the  contract  and  also  show  where  the  present  and  future  program¬ 
ming  fits  into  this  general  schemes. 
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SECTION  1 


Review  of  Vector  Spherical  Ha  rmcnics 


The  Vector  Spherical  Harmor.i^P,  T^  (0  <S>  )  form  a  complete 

J  Li  'I 

orthonormal  set  of  vector  functions  in  0  and  (j?  .  They  are  defined  in 
terms  of  the  scalar  spherical  harmonics  (  B ,  )  in  the  following 

way 
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where  the  functions  Y  (S/  ^  )  are  defined  by 
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They  form  a  complete  set  of  scalar  functions  in  Q  and^  ,  and  satisfy  the 

orthonormality  condition 
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The  functions  P  (CoS®  )  are  associated  Lengendre  functions  defined  by 
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The  vectors  e  are  defined  by 
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The  symbols  ^LM  L'M'  |  L"M"^ 

are  Clebsch-Gordan  coefficients  and  satisfy  the  conditions 
M+M'  =  M"  ^  L-L'  |^L"  LrfL' 
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The  orthonormality  of  the  TJL's  is  express 


ed  in  the  following  equation 
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An  important  equation  for  applications  is 
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The  vector  spherical  harmonics  have  particularly  convenient  differential 
expressions.  One  can  show<2>  that  the  ordinary  differential  operations 
of  vector  calculus^namely, gradient  divergence  and  curl  can  be  expressed 

in  the  following  form  _  M 

vjf (A-.-t) 
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One  can  form  another  set  of  vector  spherical  harmonics  whi  ch  we  call 

M, 
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al(  )  J  bl  (9  )  and  CL(  d  <p  )  defined  by 
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where 
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are  as  defined  earlier,  and 
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These  functions  are  orthonormal  in  the  function  sense  and  also  in 

w  M  .  M 

the  vector  sense,  i.  e.  A  L  points  in  the  radial  direction  BL  ana 

are  tangential  to  a  sphere  and  orthogonal  to  each  other  in  the  vector  sense 
for  M=M'  and  L=L'.  It  is  much  easier  to  visualize  this  set  of. functions 

•  rr  M  _  . 

in  physical  space  than  the  original  functions  \  '• 

Non-linear  terms  can  be  dealt  with  through  the  application  of  the 


following  equations 
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The  details  of. the  treatment  of  the  non-linear  terms  will  be  shown 


later. 


1.  Edmonds,  A.  R.,  Angular  Momentum  In  Quantum  Mechanics, 
Princeton  University  Press,  1957,  p.83  • 


2.  Loc.  Cit.  p.  84. 


SECTION  2 


Weather  Equations  In  Vector  Spherical  Harmonics 

I.  Basic  Equations  for  Weather  Predictions 


The  basic  equations  for  a  primitive  equation  model  summarized  below: 
1.  Equation  of  Motion: 


dc 
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(  V-V  )V  t  2Jlx  1/ 


•h  ?cP  +  JL  VP 
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2.  Equation  of  Continuity: 

+  V .  ( P V J  -  fty7  +  VP-V+  fV'V  =  o 
be'  1  *<= 


3.  Thermodynamic  Equation 

1  djj  ==  jL  i£  +  L  if> 

cf>T  dt  T  dt  p  dt 

4.  Equation  of  state  for  an  ideal  gas: 

r  -  Sr 


(2) 


(3) 
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The  variables  appearing  in  these  equations  are: 


v  =  vector  velocity  of  the  atmosphere 

p  =  pressure 

p  =  density 

T  =  absolute  temperature 

The  following  are  given  functions 


gravitational  potential 
Heat  input  to  the  atmosphere 


The  following  are  constants: 

r  =  'k  ,,  Si  =  angular  velocity  of  the  earth  =  7,  29x1  o"5  / Sec 

R  =  universal  gas  constant 

p  =  specific  heat  at  constant  pressure 

k  =  R/Cp 

=  dynamic  viscosity 

The  frict’on  term  h  V  V  is  one  of  several  possibilities  and  needs 


to  be  modified  close  to  the  earth's  surface. 


II.  Method  of  Solution  Using  Vector  Spherical  Harmonics 

The  method  used  is  to  expand  the  scalar  fields  as  series  in  scalar 

■  t«  . 

spherical  harmonics  and  the  vector  field  v  in  vector  spherical  harmonics. 


Either  the  set 


or  the  set 


AM,  BM,  may  be  used.  The  latter  provides  a  more 

L  L  L 

obvious  physical  interpretation,  so  we  will  expand  v  in  terms  of  that  set. 

The  know  fields  q,  ^  are  also  expanded  in  terms  of  scalar  spherical 
harmonics.  It  should  be  noted  that  the  expansion  coefficients  are  known. 


Ila'.  The  Equation  of  Continuity 


We  let 


1  fL  (9,<p) 

t4f  L 

2  cCL  (f,t) A* + 

h‘L  +  c”(r>t)  Cl ( <S>  cf) 


We  write: 
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We  thus  have 
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These  terms  are  attacked  in  the  following  general  n>anner: 


M  M 

The  products  A  ,  A  , 
L  L 


M  M' 

Y  ,  ,  etc.  ,  may  be  written  as 


follows:  (The  jC's  are  numerical  coefficients) 
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The  nonzero 


terms  in  VP'  ^  are  then 


I  \S  ^ 

\r  (L  ^  P**\p)  '  /  n  (20) 
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u'  .  .  yM"  (21) 
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%  p"  tL  -K,  sl’  =  2.  pi L  L 
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^  g*  •  ^  c*  =  cl  csc^L>'^r^ 
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so  that 
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The  term  '  '  f  (  V  V)  is  handled  by  the  summation  rule  for  products  of 


::)]rL„ 
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yM  The  terms  are  of  the  form 
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Tnus  the  continuity  equation  may  be  written 


The  coefficients  of  each  are  collected  and  equated  on  each  side  of 

the  above  equation  to  yield  equations  in  t  only,  of  the  form 


X  4j 

p 1  a*,  Caa  i- 
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f 
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r 

When  this  is  performed  for  all  values  of  M  and  L  we  are  left  with  a  system  of 

partial  differential  equations  involving  the  expansion  coefficients  /)  ^  , 

M  M  M 

ajj»  anc*  .  These  must  be  solved  in  conjunction  with  the  equation 

of  motion  in  the  next-section,  and  the  thermodynamic  equation  following  that. 
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lib.  Equation  of  Motion 


V  f  (7?)  7  +  iS'V  vf  +  1 
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As  b  efore,  let 
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and  given 
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f  -  i  <? :  rt 
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1 .  e  .  the 


A 

k 


M 


are  known 


is  the  z-direction  through 
earth's  axis  of  rotation. 


M, 

All  coefficients  are  functions  of  (r,  t)  only;  the  A,  B,  C,  and  s  are 
fur.rtions  of  (fy  <^)  only.  We  proceed  to  calculate  each  term  in  the  equation  of 
motion.  It  will  be  observed  that  each  terms  will  produce  one  or  more  terms  of 
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the  form 


,  when  substitution 


,  ,  ..  .  M  „  .  . .  _M  ,  .  , .  _M 

£A<r,tl  Al  +  fB(r.t)  Bl  +  fc(r,t) 

is  complete,  we  can  collect  terms  in  A,  B,  C  to  obtain  equations  in  the 

f  ,  f_,,  f  which  will  contain  only  the  expansion  coefficients  of  the 
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variables  and  known  values. 
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The  expressions  for  K  x  A  ,  B  ,  and  C 
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The  numerical  coefficients  are  summarized  more  conveniently  as 


follows : 

A 

K  x 

am 

al 
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M 

Upon  multiplying  the  first  of  the  above  equations  by  a  ,  the  second  by 

L 

M  M 

b  ,  the  third  by  c  ,  and  collecting  terms,  we  obtain 
L  L 
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From  this  it  may  be  concluded  that  the  Coriolis  term  7.SL  k  V  produces 
terms  containing 


Ahl 
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/I H 
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out  of 

each 
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rm 
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term  in 

Next  we  consider  (  V  .  V  )  v.  This  is  a  very  cumbersome  expression. 
M 

In  terms  of  the  Tj^  we  may  write  the  following  expression: 
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The  L/f^  are  numerical  coefficients.  We  may  now  substitute  for  the  TJ  m 
M  M  .  M 

-jj  »  ana  to  obtain  the  corresponding  result  for 


terms  of  A  ,  B 

L. 

the  A,  B,  C  frame.  Since  we  obtain  forms  of  the  type  aM,' D  bM"  wher 

h*  Xj  *  * 

D  is  some  differential  operator  in  r,  we  have  a  combinational  notational 
problem  which  we  circumvent  by  using  the  notation: 
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The  term  (  />  V)  /  is  then 
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where  (jf  (L',  M',  L",  M",  L,  M)  is  a  numerical  coefficient  and  D  is  some 
differential  operator.  The  principal  result  is  that  the  nonlinear  term  has  been 


M 


B 


M 


M 


C"  ’  terms,  and  the  coefficients 

Lj 


decomposed  into  a  sum  of  A  ,  , 

Ju  Li 

multiplying  these  terms  involve  only  functions  of  r  and  derivatives  of  functions 
with  respect  to  r. 


The  term  $  V /'^t  becomes  simply 


=  T  Ml  +  del  c« 

h  L  Jt  L  Jt  L  jr  L 


(41) 


-11- 


Finally,  we  need  the 


7'  c 


All  the  terms  in  the  equation  of  motion  have  been  worked  out.  However,  we  must 


evaluate 


which  is  given  by 


( 
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(47) 


-12- 


We  use  the  expressions 


r  v  r:=  z  c  *«> 

C  4  '  - 


where  o/,^  *  (M(  (ijty  4  ",  M",  ,  Jj  L  ) 


is  a  numerical  coefficient 


47  K'  =  2  C  *  cMn  C *  («) 


4-  =  2  cMei  4*  +  <Mtl  8*  .+  aJc,  C 7 


which  again  decompose  the  products  into  sums.  The  general  form  of  the 


expressions  obtained  is  therefore 


Z.  [  * V"  4  fV  (jJa, 


4>*' 
4;  M*, 
V,  J 


Similar  expressions  in  and  CM  7 

I.  J  J  J 


o<  Vf 
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The  final  step  is  to  collect  coefficients  of  the  and 

L  L  L 

and  equate  to  those  on  the  right  side  of  the  equation  of  motion.  (If  the  friction 
is  neglected  the  collection  of  terms  is  equated  to  zero.  )  When  we  do  this,  we 
obtain  for  terms  the  expression  following,  where  the  nonlinear  terms 

(  (.  V'  V)*  j  ex'  Vp  )  are  sketched  rather  than  explicitly  enumerated: 


M 

From  the  A ^  term  we  obtain,  for  a  particular  M  and  L: 

^  euK:  o  (Kci  c„  r;  ^  zsi[c,d  <^L) 

dt  Ari/'t.'  h"  l» 

*  J  J  J u 

+  CL-l  0*c  («;<■)]  +  '2l{0<",  D6  Pl'  j 


T  - i  {fa’(yw +  *  ] 


The  first  term  in  the  sum  is  the  contribution  of  &  V/fc  ,  the  second  sum 
results  from  (  ]?.  7)  /  ,  the  bracket  multiplied  by  251  i8  the 

Coriolis  term,  the  following  term  is  from  U  V  P,  the  last  term  on  the  left 
is  from  (f>  ,  while  the  right  side  of  the  equation  is  the  friction  term. 
This  yields  a  nonlinear  partial  differential  equation  with  the  angular  part 
separated  out-,  invilving  the  expansion  coefficients,  a,  b,  c,  p,  0&. 
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A  similar  performance  for  the  and  coefficients  yields  additional 

equations  in  the  expansion  coefficients.  When  the  process  is  repeated  for 
all  the  values  of  M  and  L,  we  obtain  a  set  of  coupled  nonlinear  differential 
equations  for  the  expansion  coefficients.  The  angular  part  is  separated^ 
i.  e.  the  differential  equations  involve  r  and  t  only. 
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He.  The  Thermodynamic  Equation 


The  equation  is 


—  di  =  L  'h  -r  £■  if. 

T  dt  r  dt  f  at 


or 

L  ,!$  =  ctr+  Krdj 

CP  dt  dt  f  dt 


we  use  the  expression 


c[df  =  c)jP  +  v»  (rep) 

dt  dt 


is  a  scalar  function)  to  obtain 


L,  dj  -  d_T  +  ~y.  vr  +  fcT  (dg+V  Vf>) 

cp  at  at  f 


To  obtain  only  products,  we  multiply  by  p  and  obtain: 

fit  =  f  flr  4-  v  •  (7T)  )  +  kr(  dj  t  v-vf) 
dt  \  dt  x  ot 


As  usual,  we  expand  the  unknown  functions  p  and  T  in  scalar  spherical 

harmonics  and  the  velocity  in  terms  of  vector  spherical  harmonics.  It  will 

be  noted  that  manipulations  for  p  (  +  v,Vt)  and  k  (  df  +  if  )  are 

^  t  Jt  i 

formally  the  same,  so  only  one  will  be  indicated.  These  terms  involve 
V  >  (  V  )  which  we  write  as 
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i /Xvr)  = 


s(ia:^fki  ^>■(’’(2  K'rpjj 

\  *—  L  U>L> 


Hl 


i  2  i,y^<'*r^:c) 

*,L 


(57) 


As  in  the  case  of  the  equation  motion,  the  nonzero  terms  that  result  are 


M  ;M'  bm'  CM'  .  BM'  .  Using  the  decomposition 

r  •  T  I  »  r  •  t  «  ’  I  ,  *  T  , 


‘L  *  "L'  ’  ~L  *  L' 

rules  for  these  terms  we  obtain: 

V.VT  = 


z«;«  n?  <&.  n:  *  i  *:  % c*  c  *  z  <?<« 


'ic 


(58) 


y/W'/ 

'  4“ 


This  term  must  now  be  multiplied  by  T  which  can  also  be  handled  by 


summation  rules.  We  have 


f 


■  (  f.vr) 


( 


i  f  r j(z 

£/ii 


z  «>rrr» ; 


(59) 


which  will  yield  terms  involving  sums  of  products  of  coefficients  of  T, 


? 


c*  t"'  *  r^'  rN“"  l" "> •••  i ^ 

6  L>  •  T  L"" 

^  IT  iLr\\ 


(60) 
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I 


also,  we  will  encounter  terms  of  the  form 


and 


.a/"'  ' 

U 

Ca  a  (  M 

l  / 

which  can  be  collected  for  specific  values  of  the  indices  M"M,  L,,M.  This 

! 

operation  will  not  be  carried  out  in  detail  because  of  the  complexity  of  the 
algebra.  1 

dT/z>t 

is  also  handled  by  the  tie  composition  rule;  we  have 


The  term 


P 


fir.. ,  (i  fMr:  )(z  £z"'  rj) 

1  dc  J  1  MX  it  ■  '  C  ' 

\ 

,  uli  L>' 

M,0  u 

*  jU  J 


(63) 


(64) 


i  i 


The  left  side  of  the  equation  is  also  handled  by  the  decomposition  rule  to  obtain 


2  fl  iv  .K- 

lA/V...  U  '  - 


(65) 


i  / 


W 


where,  the  q  are  given  coefficients  repres enting  heat  input. 
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I 


Finally,  the  term  kT 


C  \S  ^  ^  f  )  is  handled 


analogously  to  the  p  ^ +  \/  .yf  j  term. 


Upon  carrying  out  that  operation  and  collecting  all  the  coefficients  .of  the 
M 

Y  's,  we  have  a  set  of  coupled,  nonlinear  partial  differential  equations 


involving  T 


M 


M 

p.  ,  and 


M 


,  and  as  usual  the  angular  part  does  not 


L.  L  L 

appear;  so  that  the  equations  are  of  first  order  in  r  and  t. 
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SECTION  3 


Pressure  Coordinates 


If  one  wants  to  use  pressure  coordinates  in  the  spherical  representation 
r  is  replaced  as  an  independent  variable  by  p.  The  following  treatment  is 
incidentally  valid  independently  of  the  specific  meaning  of  p,  it  can  be  any 
scalar  expressed  in  spherical  coordinates.  Let  us  apply  equation  (35)  in 


the  section  1.  to  £(0^^  where  p=p(r,  <p  ). 

and  similarly  for  y,  and  z. 

a  x '  “*jp  ^ 


Then 


Thus, 


||af 


and  (35)  becomes 


+  ^\lzr7S~  MT  I  -  H  l  m>  f<>) 


(1)  can  be  used  to  find  the  appropriate  differential  equations  in  pressure- 
coordinates,  but  sometimes  there  are  more  direct  ways.  The  gradient, 
divergence  and  curl  expressions  become 


20 


The  continuity  equation 

can  be  handled  as  follows. 


P 


V ■  V+  V. 


vr 


=  o 


Substitute 


_  I  dip 

f  ^  lnto  the  continuity  equation  to  obtain 

4!? ■*  v:?|£  =  0 

5  K  -  v.TPi,|rr.v*-^.Tl,=  o 


Write 


z>t 


2/z. 

L*-)  and  using  Eqs  3,  4,  5 


3uo 

7>K  ~r  2>/L  |_ 


7.vt 


af 


.  2?.7p 
r 


but 


*t\Y '2£?*.A%y*}’Lvt.b\y  $s?  vk  c*"  *  ^  vl 

bn\r*f>  1  L  T  - ^ /-=ajT  T' 


Therefore, 


<76t> 

2>A.  &/2. 


&zuk<f>-K<p>iy; 


S3 
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and 


then 


to  M 
if 


(9) 


where  the  -^..’2  - 

1>L 


pafF 


substitution  has  been  used. 


Note  that  we  have  kept  the  vertical  velocity- component  in  the  derivation. 
If  we  leave  it  out  the  continuity  equation  reduces  to 


-  %  iT(p) 

a  p  ~  u  1 


(  10 


Thus  one  can  see  that  only  part  of  the  horizontal  wind  contributes  to 

namely  that  part  in  the  expansion  which  depends  on  functions  (©,<p) . 

The  zonal  wind  which  depends  on  terms  of  the  form  C>,f)  and  Rossby 

M 

waves  which  depend  on  terms  of  the  form  CL  )  do  not  contribute  to 

.  This  can  perhaps  be  an  alternative  explanation  to  the 
smallness  of  (jJ,  namely,  the  large  features  of  the  atmospheric  circulation  do 
not  contribute  too). 

N  ote  that  equation  (1)  shows  that  there  are  no  derivatives  with  respect 
to  r  left  when  p,  ©  and  are  the  independent  variables  so  in  this 
formulation  it  is  not  necessary  to  use  the  hydrostatic  balance  condition  to 

make  the  substitution  \  \ 

=  -  0  C*. 

7>K.  (  J 

£3 
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Note  also  that  r  still  remains  in  the  equations  and  it  does  not  seem  possible 
to  eliminate  this  variable  from  the  equations  which  means  that  one  cannot  really 
cast  the  equations  in  pressure  coordinates  at  all.  This  does  not,  however, 
appear  to  be  too'  serious  since  putting  r  equal  to  the  radius  of  the  earth  is 
a  good  approximation. 


SECTION  4 


The  Transform  Method 


In  recent  years  a  transform  method  (1,  2)  has  been  developed  which 
appears  to  make  spectral  methods  competitive  with  the  finite  difference 
method  as  far  as  speed  is  concerned.  We  give  a  brief  discussion  of  the 
basic  ideas  and  show  that  they  also  work  for  vector  spherical  harmonics. 

Consider  a  differential  equation  for  a  scalar  field  Ip  0  expre 
in  terms  of  the  independent  variables  t,  Q  and  ^  .  Let  Ip  satisfy  a 
partial  differential  equation  of  the  form 


ssed 


+  L*  =  N(f) 

£  U 


(  1  ) 


When  L  is  a  linear  operator  ana  N  a  non-linear  term.  L  and  N  may  or  may  not 
contain  differential  operations.  Let  us  ignore  the  linear  term  since  it  is  easy 
to  deal  with  and  just  write 


-  N  ( O')  <  2 1 

at 


The  conventional  spectral  method  appro? ch  to  solving  (2)  is  to  expand  <p  in  terms 
of  a  complete  orthonormal  set  of  functions  in  0  and  ^  ie.  scalar  spherical 

harmonics  YT 

Vi 


(  3) 


ML 
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To  be  specific  let  N  ( 
dimensional  reasons. 


)  be  equal  to  9C  y  where  a 
Thus  (2)  becomes 


is  a  constant  used  for 


dt 


(4) 


and  substituting  (3)  into  (4)  we  get 


Now  product  term  of  the  form 


can  be  expressed  in  terms  of 


sums  over  individual  terms 


y 


by  the  use  of  the  formula 


^■yo]  <  Lolo|iTo^<lVlV|lm+m|>  YjT 


Let  us  abbreviate  the  result  of  substituting  (6)  into  (5)  by  expressing  (5)  as 


(7) 
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£6 


's  we  multiply  both  sides  by 


New  using  the  orthonormatity  of  the  V’ 

u 

LYcT s in  ©  d  ed <p 


M 


and  integrate  Then: 


p  r  ^  m"l1'  c  L  c:*  m+m'Ip  L,i 


(8) 


or  • 


d  CL, 

dt 


*  =  a-  2]  a“  a£‘Mi (LY  ^.Kp- M') 


m'(l'  iT 

i  i 


(9) 


The  basic  idea  is  to  replace  the  time-consuming  calculation  of  the  right- 
hand  side  of  (9)  by  a  faster  method. 

Let  us  rewrite  (7)  in  the  form 


mW’ 


(10) 
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1  ' 


Z7 


Multiply  both  sides  by  [V>,<pA  S\W0o(0cf  ^  and  integrate  (10)  to  obtain 


— ^ 

dt  K 

where  the  explicit  form,  which  one  does  not  want  to  use,  for  the  right-hand 

side  is  that  of  the  right-hand  side  of  (9).  Instead,  there  is  a  less  time 

consuming  way  of  finding  d  and  to  attack  (10)  in  the  following  way. 

i  n  in 

Suppose  we  truncate  L,L  ,L  and  L  at  N.  Abbreviate 


,  mW 


D  £  twc,X«,  =  Jjw) 

ll"s0  M  M  +  M<L  L 


Write 


YjV.p)  =  Xje)  e 


then  form  (5) 


* *  fctfcWfc'WM 

=  v  j 
.  mwjl;'  is 
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Now  let  6 
Let 


and  take  on 


discrete  values 


and  insert  them  into  (14) 


■  <P  =  2.TT  1C 
'k  4  M 


Multiply  each  side  of  (14)  by 

-2JTL  rvxk/ 

__1_  e  /4  n 

4N 


and  sum  over  k 


Z,  _2TT  i  rvi  l<  / 

/4'1 


d  X  (e.)e  4N 

m'+m",  ii"  r  i 


2TTrv» 


wk’ 
N  . 


(15) 


Now 


4-N 


-  7 

4-N  L_j 


2. it  £,  (m+m'L 


-  5“ 


/<=  * 


M  *  m 


(16) 


X%)sg(») 


-5. 
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m  is  fixed  so  the  sum  over  m  vanishes  and  we  have  a  set  of  equations,  one 
tor  each  m#  (17)  then  become 


£  -  w 


■(18) 


Xjj"  (  )  can  then  be  considered  a  set  of  matrices  labeled  by  m  with  Li  and  j. 

as  the  indices.  Finding  the  inverses  ci  the  set  of  matrices  X  we  can  solve 
(18)  for  d  which  is  what  the  problem  to  be  solved  was. 

We  will  now  indicate  how  the  transform  method  is  applied  to  vector  spherical 
harmonics.  Since  we  have  an  analytic  e-'oression  for  the  term  (  v  •  ^)"^in 
the  vector  spherical  harmonics  formulation  and  will  be  able  to  use  it  to  check  the 
accuracy  of  the  numerical  computation,  we  use  (V«V)  v  as  an  example. 

The  key  to  putting  the  transform  technique  to  practice  is  to  make  sure  that: 
all  V  operations  are  performed  in  coefficient  space,  all  non-linear 
operations  are  performed  in  physical  space.  Of  course,  one  must  make  sure 
that  the  density  of  physical  space  grid  is  consistent  with  the  number  of 
coefficients  and  the  particular  non-linear  operations  being  performed. 

Consider,  thus,  the  nonlinearity: 

(7.V)7 

It  can  be  expressed  as 

( 7.  V)  v*  =  i  V(v’.v')-  v*x.(Vx  7) 


Remembering  that  we  will  be  computing  the  vector  spherical  harmonic 
coefficients,  v  ,  of  the  variable  v  and  that  in  coefficient  space,  the 
curl  of  an  expanded  variable  is  easy  to  form,  vx(  V  X  V  )  can  easily  be 
formed.  It  requires: 

It  requires: 

o  forming  curl  in  coefficient  space  and  transforming  to  physical 
space.  (  7x7  )., 


m 


o  transforming  the  v  to  physical  space.  (  V  ). 


o  forming  the  cross  product  at  each  point  on  the  physical 
grid  (  7x  (V  x7))  . 

The  gradient  of  v.  v  is  a  little  more  difficult.  In  the  space  defined  by  the 

c  •  unit  vectors,  v  can  be  expressed  as: 

n 

7?  A  A 

v  =  v  e_  +  v.eo  4  v+e+  - 

and  using  the  dot  product  rules  in  the  e  space 

V.7:  -2VV  +V.1 

-  4  O. 

and  then  applying  the  gradient  operator  produces 

£7(7.3)*  v  7vo  -  v  -  v,Vv_ 

It  can  be  shown  that  for  a  real  vector  "v^ 

v*  =  _  V* 
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I 


1 


leading  .to  the  relationship 


Y  Vv_  = 


(v.VvV) 


The  remaining  steps  in  computing  the  advection  term  accomplish  an 

i  1 

evaluation  of 

-1  vVv  -2  Re(vVv+)  . 


followed  by  a  VSH  expansion  of  the  non-linearity. 


o  forming  gradient  of  the  e^  component  of  v  and  transforming 
to  physical  space.  (  W  ). 


,  J  ,  t  A 

o  forming  gradient  of  the  component  of  v  and  transforming  to 
physical  space.  (  V  V+  ).  ,  i 


r 

o  forming  the  product  of!the  above  gradients1  and  the  appropriate 


components  of  v  over  the  physical  grid.  (  V  V  V  and  vVv  ).' 


a  o 


o  Sumining  the  above  results  over  the  physical  grid  forming  ((vl^)"v). 


°  ,  Transforming  the  advection  term  to  coefficient  space  producing 

,  i 

the  VSH  coefficients  of  (  • 


■  8- 


3; 


.  .  . . . . . . 


To  illustrate  the  advantage  of  using  the  transform  technique,  consider 


the  vector  example 


where 


— *  v  M  — p-  M 

v  “  iVT:i 

M.L 


Formal  substitution  produces. 


<3v  M 

_ -  3| 

at 


LaU  '  T".*l  j  L  J  L 


lA  3  L  T  L 


Assumming  that  the  C's  are  precomputed  and  available,  a  rough  operation 

count  (where  L  is  the  truncation  value  of  L)  is: 
m 


+  C.V.V  —  6  real  mpy  +  4  real  add  (v  complex,  C  real) 

—  5  operations  (an  operation  is  mpy  +  add) 

M  M 1  a  i  Z 

the  number  of  V  =  the  number  of  V  t  ,  — 

jl  Yi:  — r  • 


the  number  of  V 


(superscript  determined  by  M&M'  above)  =  3  L 


»  »  total  number  of  accumulations  =  total  number  of  C's  to  be  stored 


and  the  total  number  of  operations  ~  135[? 

4 

Approximately  80%  of  the  C's  included  in  the  computation  above  are  equal 
to  zero.  Assumming  that  these  could  be  bypassed  at  no  cost,  the  totals  become 

number  of  C's  to  store  ~  1.35  IT 


and  number  of  operations 


~  6.75  t! 
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Use  of  the  transform  method  involves  the  following  steps: 


1) 

2) 

3) 

4) 

5) 

6) 


Compute 

Compute 

Compute 


Over  the  physical  grid. 
Essentially  four  vector 


Compute  V 

Combine  over  physical 
Transform  (7.  V  )  V 
C£  5  L.  operations. 


transformations . 

operations  each 

rr'  2 

grid  to  form  (  V.  operations 

to  coefficient  space  (one  vector  transformation) 


The  total  number  of  operations  by  this  technique  is  approximately  25L 
This  represents  a  significant  savings  for  large  (see  Table  below) 
There  is  also  a  large  savings  in  storage  locations. 


Operation  Count  for  Computing  the 
'  (  V.V)  V  Portion  of 

5  4  3 

Straight  Forward  Method  .  .  .  .  6L  m+  12L  m+  6L  m 
(assumes  existence  of  selection  rules  that 
eliminate  80%  of  the  terms  at  no  cost) 

3  x  2 

Transform  Method  ....  25L  m+  (415  +  200  log^L^)  1^  +  251^ 


L 

4 

8 

16 

32 

6.4  .  . 

128 

m 

SFM 

QM 

.65 

3.  2 

17 

94 

500 

2500 

T~f 
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SECTION  5 


The  Mintz-Arakawa  Model  in  Vector  Spherical  Harmonics 

A  simplified  version  of  the  Mintz-Arakawa  model  has  been  written  in 
terms  of  the  coefficients  of  the  vector  spherical  harmonic  series,  representing 
the  dependent  variable  T.  Also  used  are  the  scalar  series  of  the  scalar 
spherical  harmonics ,  used  for  temperature,  ground  pressure,  etc. 

Left  out  in  this  preliminary  version  in  agreement  with  the  statement  of 
work  in  the  contract  are  moisture  and  radiation. 

In  order  to  accommodate  the  continuous  nature  of  the  model,  the  friction 
term  at  the  bottom  has  been  modified.  Furthermore,  it  is  assumed  that  the 
earth's  surface  temperature  is  given  as  a  function  of  position  and  time.  The 
surface  temperature  can  lateron  be  calculated  from  zero  heat  balance  on  land, 
if  desired. 

The  Mintz-Arakawa  Model 

The  Mintz-Arakawa  equations,  as  given  in  [3]  are  the  equation  of 
horizontal  motion  (  coordinates) 

J^(ir  v)+  (V.7r7^)7  )  +  f  T  X 

-h  (TTrtxivTr)  =  tt*  F 

The  thermodynamic  equation 

-  T) +  *  (1rcrTr)  T&) 

—  7TCX  +*  G*  ^  VTT-h  7T^)s  7T  // 
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3G 


The  continuity  equation 


^  t 


The  moisture  continuity  equation 


(-rr%)  +  v.(rrf.  'r<? 


7?7~ 


V 

The  equation  of  state 

/°~P 

and  the  hydrostatic  balance  equation 

-f  7T  OC  ar  a 

<r 

In  these  equations 

=  surface  pressure--top  pressure 
=  horizontal  velocity  (in  (T  coordinates) 
where  p  is  the  pressure 

A'Ar 


TT  =  Ps-PT 


o' 

C T 


k 

f 

-Q. 

* 


d£* 

dt 


=  vertical  unit  vector 
_  2-0.  sin  JP  where 

=  rate  of  rotation  of  the  earth  and 
=  latitude,  positive  worthward  from  the  equation 

-  geopotential 

/ 

=  specific  volume  =  j£~ 

-  horizontal  vector  frictional  force  per  unit  mass. 
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specific  heat  (for  dry  air)  at  constant  pressure 


T  =  temperature 

• 

H  =  diabatic  heating  rate  per  unit  mass 

q  =  mixing  ratio 

^  =  longitude 

=  rate  of  moisture  addition 
The  boundary  conditions  are: 

At  /  over  land 

cr  =o 

*  -  geopotential  of  the  earth's  surface 

and  F  =  0  where  F  is  the  vertical  heat  flux  at  the  surface. 

H  H 

At  CT  =1  over  ocean 


i>  - 


Ts(M) 


and  at  p  =  pT  : 

Sit  =  o 
cit 

which  is  the  fr  surface  condition. 

These  equations  are  written  in  momentum  form>that  is,  the  variables  V,  T  and  q 
are  multiplied  by  IT  ,  so  instead  of  V,  T  and  q,  the  variables  are  TT 7,  TTC  Tj 


This  can  be  accomplished  by  multiplying  the  equation  of  motion  with  Tf*  , 
the  continuity  equation  by  v  and  adding  the  two  resulting  equations.  In  the 
next  section  a  short  derivation  of  these  equations  is  given  and  the  major 
approximations  made,  are  listed.  However,  the  equations  will  be  written  in 
the  original  variables  v,  T,  etc.  rather  than  77V  etc.  The  reason  for  this  is 
that  besides  the  variable  (  TT  v)  also  the  quantity  v  is  used  (see  for  instance  the 
equation  of  motion  where  the  second  term  on  the  left  is  7*^7V  J  V  .  This 


quantity  v  is  obtained  from 


irv  ) 


,  which  requires  a  division.  It  was  thought 


to  be  desirable  to  avoid  this  at  this  point.  Also  the  reasons  for  using  etc. 
are  mainly  that  in  the  finite  difference  approximation,  it  is  necessary  that 
momentum  is  conserved  by  the  mathematical  approximation  techniques  used. 

This  simplies,  in  that  case,  to  use  (  TT*  V  )  instead  of  (/"  . 

The  following  section  reviews  the  derivation  and  lists  the  major  approximation 
used.  In  that  section  the  equation  for  the  moisture  has  not  been  reviewed  as 
moisture  has  been  deleted  from  the  Vector  Spherical  Harmonics  Model  at  the 
present  time. 

The  following  approximations  are  made  to  obtain  the  quasi-static 
approximation.  [1,  p  20-22]. 

Subscript  h  indicates  horizontal  component 

1°  Vertical  acceleration  is  neglected 

2°  Horizontal  component  of  ./l.  is  ignored  in  the  expression  for 
the  Coriolis  force 

3°  (f/ih  i”  by 


3  (57ih  is«pla«dby  ("5“ 

4°  div^  V  )  i®  replaced  by  div^  ^0 
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The  equation  of  motion  can  be  written  then  in  two  parts.  If  W  is  the  vertical 
velocity  and  the  subscript  h  indicate  the  horizontal  component: 


/  }  'Vj  —  —  )  h  .  /  —  \ 

( Friction- 


( i ) 


which  is  the  horizontal  part.  The  vertical  part  becomes 


ht—v 


(  2  ) 


The  friction  term  of  the  horizontal  part  will  be  specified  later;  in  the  vertical 
part  it  has  been  neglected  [l,  p  33-41]. 


Furthermore,  the  continuity  equation  is  written  without  approximation: 


while  the  thermodynamic  equation  becomes: 


4#  +  7  B  4-  w  IS-m  yr  (&■)*<? 

where  0^,  is  the  potential  temperature  ^  /°°.  J  ^ 

<2-  -  C..  P 


^  C-M- 

k  =  £  i. 


U  ) 
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Pressure  Coordinat^n 


in  these  coordinates  pressure  instead  of  r  is  used  as  the  vertica,  coordinate. 
This  simplifies  the  equations,  though  it  complicates  the  boundary  condition. 
However,  it^also  agrees  with  customary  meteorological  terminology  -  Any 
function  Y  <*.  y.  e)  will  be  changed  to  a  function  ^  <x,y,  p(x,y>  a|- 


if  co  =»  a 

c/t  bf 


+-  .  V /* 


the  continuity  equation  ca»be  given  in  the  new  coordinates  as 

O  GJ  /V  &  a  J  „ 

This  is  the  continuity  equation  in  pressure  coordinates. 
Equations  (1)  and  (4)  become 


l  5  ) 


Friction 


rf*  i  ?4  * 

Z6  / 


Furthermore: 


where 


2i  _  e<  _  ^r 


*  =  2Z 
/> 

-6- 


=  heating 


(  6  ) 


(  7  ) 


(  8) 
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C*  Coordinates  [i] 


In  order  to  accommodate  Orographic  features  and,  in  general,  make  the 
earth's  surface  a  surface  on  which  only  )  vary  a  new  variable  is 

defined: 


(Tag  JL. 

TT 


(  9  ) 


where  is  the  surface  pressure.  Also  from  here  on,  V  indicates 

the  horizontal  component  of  the  velocity,  unless  otherwise  indicated.  The 
following  are  the  transformation  rules: 


(h) 


—  ^ r  =  s  2— 
2c r  "  7T  )cr 

(U  '  ’ 


Now  /  )  __ 

1  > 7r}tT 


*<r  ' 

v/t. 


d  /T>(T  \ 
■>C T  ^ 


—  —  which  follows  from  eq.  (  9  ) 


Tar 

~>x 


(  10  ) 


Then  in  vector  form: 


(  H) 


i_  /2r )  In 

14 

T_  22ZL  ^ 

’7 r  I*  TP 


(  12) 


'■•'""'iiiiijjai 
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First  the  continuity  equation  will  be  changed  to  the  new  coordinates. 


W  e  had  the  continuity  equation: 


The  subscript  C  will  be  left  out,  from  now  on  so  that  becomes  cu 
Using  eqs  (  10  ),  (  11  ),  (  ^  ),  equation  (  13  )  becomes: 


^  (T7  Cj) 


~UL+.rr‘^-(lZU(r71 


2><r 


V  IT  -h(TV7r. 


2  \/ 


4-  1/.  7  7T+  TTV.  ^ 


sO 


as 


iLZ?* ——  O  because  ^is  not  a  function  of  G~  .  (Earth's 


surface  is  0“  =  1  surface). 

So  the  new  continuity  equation  is: 


2JZT-I-  v./'m/)  -t-2£Z^  =  o 
h  ■£  %  (T 

Using  eq.  (  /o  )  the  equation  of  motion  becomes: 


(  14  ) 


FRICTION 


(  15  ) 


f 
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with  the  hydrostatic  equation: 


2A  =  _  ZT 

2  <r  C~ 


(  16  ) 


Integrating  eq  (  14  )  from  =  0  to  (T  -  1  gives: 

^21  —  V.  f(z rv)cL<r 

4  ' 


where  use  is  made  of6J  =  0  at  6r=  0  and  0*=  1 

75T'. 

Eliminating  —  in  eqs  (14  )  and  (  17  )  results  in: 

&T 


(rrv)+j  v.(irV)  trj 


(  17  ) 


(  18  ) 


Noting  that  in  equation  (  14  )  0  gives 


"2>  co 

!><r 


=  -£[} T  +  7'CTrv)] 


which  is  equivalent  to  (  14  ) 


(19  ) 


Eq  (18)  can  be  used  to  calculate  OJ. 

The  thermodynamic  equation  becomes: 

V  +  ^  v  e  +^^=/r/-4-jV 

h{  7  r  2  <r  <£,  1 trir'i 


(  20  ) 
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Equations  (15),  (17),  (18),  (20)  jare  to  be  solved  as  follows  [2,  p  341-344] 

The  initial  data  to  be  given  are: 

Tf(&,  f,  ‘=0) 

T(b,  t=oi 

vr(6'  f>  t=0> 

First  from  eq.  (  16  )  calculate  .  To  this  purpose  must  be  given. 

This  is  orographic  data.  Next  W is  calculated  from  eq  (is  ).  This  is  the 

diagnostic  part  of  the  computation. 

— 

N ext  with  the  help  of  eq  {  15  ),  is  computed.  Eq.  (17  )  is  used  to 
calculate  4“^.  Eq.  (  20  )  gives  £j§r.  The  results  are  used  to  calculate 

_  ot  ■  bt 

the  Vj  TT* »'  T  at  time  t=t  +  L  t 
Friction  Term 

For  a  tentative  description  of  the  vertical  eddy  transport  the  expression,  given 


by  Robert  [5]  was  used: 


Friction  =  T7  (fi 


2  f  h>z  ^  ^ 

TjT 


0(.  =,  0002  hour 


Boundary  Conditions 


If  eddy  viscosity  and  bottom  friction  are  neglected,  sufficient  boundary 
conditions  are 


CJ  =  0  at  0“  =  1 

CU  =  0  at  (T  =  0 
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The  first  condition  specifies  a  zero  mass  flux  into  the  earth.  The  second 


dp 

condition  is  interpreted  as  the  free  surface  condition  ” 


0.  Though  the  i 


condition  at  the  top  and  bottom  are  the  same,  they  mean  different  things, 

namely,  at  the  top  ^  =  0.  Now  p  =  & Tf  so  rj“  =  *jJ7  +  ^  As 

(7"  =  0  at  the  top  and  it  follows  that  0J ,=  0  at  the  top.  i 

dt  i 

i 

These  two  boundary  conditions  are  still  sufficient  if  only  the  eddy  viscosity  is 

‘  i  |  '  I  i 

neglected.  This  is  because  of  0J=0  at  the;  two  boundaries  and  the  occurance  of 

only  in, the  equation  of  motion,  s.  g.  (15).  [2,  p.  96,  97] 

In  those  cases  V  the  surface  velbcities,  could  be  calculated  from  the 
S 

differential  equation  (  15  )  by  putting  OJ  - —  =  0  and  using  a  suitable  expression 

.•  a<r  iii 

for  the  surface  friction.  In  many  cases  this  frictional  force  is  assumed  to  be 

I 

very  large  so  that  is  effectively  equal  to  zero.  If  this  is  not  the  case 

i  i  i  i  •  i 

must  be  calculated.  To  do  this  with  the  differential  equation  might  be  top 


inaccurate.  The  forces  that  are  acting  on  the  air  particles  are  dominated  by 

•  '  i 

frictional,  viscous  forces  which  are  not  included  in  eq. 1  (15).  These,  forpes 

'  1  1 

work  in  a  layer  of  about  1  mm  thick  [1,  p  38],  ,  Within  this  layer  the  wind 
velocity  changes  rapidly  to  become  zero  on  the  ground.  At  the  top  of  this 

i  , 

layer  V  is  no  longer  zero.  It  is  therefore  better  to  give  some  boundary 
condition  on  v  in  that  case.  The  same  is  true  if  the  viscosity  is  no  longer 

l 

neglected.  The  vertical  transport  of  eddies1  can  be  taken  into  account  with  a 
friction  term  K  S — £  [1,  p  38,  39]  in  which  case  the  equation  becomes  2nd 

o  X . 

order  inZ  and  two  boundary  conditions  liave  to  be  given  [2,  p  96]  Rather  tl^an 


i 


putting  the  horizontal  velocity  equal  to  zero  at  the  surface,  a  drag  coefficient 


is  introduced  [5], 


(1*)  =-f'7 
^  h  cr  '  $  s 


where  <£  is  an  experimental  constant.  Robert  gives  £  =  40  which  corresponds 
to  a  surface  layer  of  about  800  meter 

At  the  top  of  the  layer  —  =  0  [2,  p  96],  These  are  the  conditions  used 

here.  However  more  work  on  this  problem  will  have  to  be  done. 

Vector  Spherical  Harmonics 


If  vertical  acceleration  is  neglected,  the  equation  of  motion  falls  into  £ 
parts,  one  part  becomes  the  equation  of  motion  (  15  )  for  the  horizontal  motion, 
the  other  part  becomes  the  hydrostatic  equation  (  .16  ).  Consequently  V can  be 
expanded  into 


v'  =  L.  \n>  ■%  o  C. 

1  ci  -  /  J 


(  21  ) 


where  the  and  are  the  vector  spherical  harmonics  as  defined  in  [4], 

.nM 

They  are  functions  of  the  coordinates  (&.P)  only.  The  and  tit 

'  4L  CL 

are  functions  of  ( ^  ).  They  lire  scalars  and,  in  general  are  complex. 

M 


Similarly  the  vertical  (T  -velocity  is  written  as: 


CO  =.  2_  CO  (<r;0  yL 


(22  ) 
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the  ground  pressure  as 


the  geopotential 

7*  =  lM 


the  temperature 


(  23  ) 


(  24  ) 


(  25  ) 


the  potential  temperature 

o  =  H  y*  < 26 

T  L,M  Tl  ,L 

Before  substituting  these  expressions  into  equations  (15  ),  (  16),  (  17  ), 

(  18  )  and  (19  )  it  should  be  pointed  out  that  the  use  of  transform  methods  avoid 
the  necessity  of  calculation  coupling  coefficients.  Vector  spherical  harmonics 
expansions  of  non-linear  expressions  such  as  TTV and  (V,  v" can  readily  be 
obtained  without  any  further  approximation  except  truncational  at  maximum  • 
value  of  L  in  equations  (21  )  to  (  26  ).  The  coefficients  of  the  VSH  expansions 
of  (/.V)  ~  for  instance  are  indicated  by 


In  other  words 


v) 

^  =  X  /? 

\vv 

7  A"- 

V 

Lh 

C" 

c£.  L 

An  exception  1 

own  symbol 

was  made  for 

Cp”9]  „ 

A  q' 

which 

ft 

was  given  its 

cy  M 

=  Z 

VA 

- 

V  v 

r 

7 

To  change  equations  (  15  ),  (16  ),  t  17  ),  (18  )  and  (19  )  into  equations 

for  the  spherical  harmonics  coefficients,  use  is  made  of  the  expressions  where  the 

vertical  velocity  has  been  ignored. 


II 

1/42*/  /J3  * 
a.  /  1 

(27  ) 

y^ci+s  ^  yfi 

(  28  ) 

s  •/<?*- 

o 

(  29  ) 

Also,  in  [4]  a  general  expression  is  given  for  >C  V  .  From  this: 

^  AV,  -  T07T  &  ^ 


X vl, - 


(30  ) 


(31  ) 
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Substituting  eqs  (  21 )  to  (  26  )  into  e,  (  15  )  and  making  use  of  eqs  (27  ) 

M  ,  thp 

to  (  32  )  and  collecting  terms  multiplying  the  same  an  L 


following  equations  replace  eq.  (15  ): 


(r'i-XL 

^  Ixr 


(  15b  ) 


The  hydrostatic  equation  becomes: 

i nr  ■  4 


(  16a) 
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Equation  (  17)  becomes: 

Itt” 

fa*  « 

1i 


while  the  equation  for  the  vertical  CT  -velocity  becomes: 

A  division  by  has  to  be  made  here,  however,  this  is  done  in  an 
expression  for  a  small  quantity. 

The  thermodynamic  equation  is  changed  to: 

*  ATZ  ms  e*2*l  +/uli/  = 

L  LM  G*  TL  l  Jyji  l 

Radial  Dependence  ' 


Hi 


(17a  ) 


(  18a  ) 


(20a  ) 


The  equations  (  15a  )  to  (  20a  )  still  contain  the  derivative  and  the 

integration^/*e/^flls  the  etc  are  not  eigenfunctions  of  the 

Navier  Stokes  equation,  no  natural  radial  functions  have  yet  been  decided  on. 
Also  it  should  be  kept  in  mind  that  the  amplitude  of  these  radial  functions  are 
functions  of  time  to  be  evaluated  by  solving  the  differential  equations.  The 
result  is  a  set  of  coupled  ordinary  differential  equations  with  these  amplitudes 
as  the  dependent  variables  and  time  as  the  independent  variable.  The  simplest 
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functions  are  the  powers  of  .  They  have  the  important  property  that 

products  of  two  powers  is  again  a  power  of  (7*  .  This  will  simplify  the 
treatment  of  the  non-linear  terms  considerably.  Therefore  power  series  will 
be  used  as  a  first  try. 


Now  the  spacial  boundary  conditions  impose  certain  conditions  on  these 
coefficients  to  be  maintained  at  all  times. 

Power  Series  Expansion  in  0** 

Suppose  we  consider  here  one  possible  approach.  Let 

A  typical  boundary  condition  would  be  v  =  0  at  the  surface,  where  ^-1.  So 
the  boundary  condition  becomes: 


V  „  r-  M(k)  _  t 

2-  V-  =r  )  _  ^  =  O  :  CTts  f 


r 

Z-  T/  = 

U  U=o  <-L 

and  the  coupled  differential  equations  have  to  be  solved  with  these  extra  condit  ons. 
The  method  adopted  will  be  described  using  power  series,  however,  it  is  not 


necessarily  limited  to  power  series. 

Suppose  emperical  evidence  suggests  an  Nth  degree  curve  for,  say,  the 
temperature.  So  there  are  N  +  1  coefficients.  If  there  is  one  boundary  condition, 
for  instance  T  =  o  at  (T=  1,  it  means  that  an  Nth  degree  curve  must  be  fitted 
with  the  condition  that  T  =  0  at  T  =  1.  This  will  give  N  degrees  of  freedom. 

If,  on  the  other  hand,  the  temperature  is  given  in  N  points,  (besides  T  -0  at 
(T=  1),  the  Nth  degree  curve  can  be  reconstructed  at  any  time  without  losing  any 
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information.  Instead  of  using  the  (N+l)  coefficients  as  variables  it  is  therefore 
equally  possible  to  use  the  value  of  T  at  N  levels  as  variables^*  the  value  of 
T  at  the  boundary  can  be  obtained  from  the  boundary  condition  and  the  N  values 
of  T.  Furthermore,  there  is  another  advantage.  Writing  out  the  equations 
for  both  methods,  the  equations  for  the  coefficients  get  extremely  complicated 
while  the  equations  for  level  values  remain  the  same,  no  matter  how  many 
levels  are  used. 

The  following  method  has  been  adopted  to  work  with  the  layer  values  of 
the  variables,  and  still  retain  the  advant age  of  continuous  (T  dependence. 
Suppose  that,  as  an  example, 


Initially  the  VM  are  calculated  from  the  initial  conditions  for  /  .  Suppose 
CL 

"™"l*  Ivl  (7* 

l /  is  given  in  N  layers.  In  that  case  V™  is  known  at  N  different  (T  's. 

If  ,  in  a  particular  model,  is  to  be  represented  by  a  Nth  degree  curve, 

CL 

in  N  +  1  layers  have  to  be  given  initially  to  compute  the  coefficients 
,  If  a  boundary  condition  is  given,  say  ^  ^  =  0  at  (T~  1,  thenj^"^/^  ss  O 

and  N  layer  values  are  sufficient  to  calculate  the  constants.  The  boundary 
condition  provides  one  more  equation.  From  eq  (  33  )  the  N  equations  are 

4r',=L<S(i><rl 

X  tL  * 


(  34  ) 


where  j  indicates  the  layer.  The  coefficients  CT*  a re  the  known  pressures. 

The.  VM(k)  ,  " 

4  L  ’  Canbe  calculated  with  the  conditiun  (in  the  example)  that. 


„  *1*1  „ 
At  ~  ° 


The  solution  can  be  written  as: 


-  z 


k* 


(  35  ) 


Which  is  a  linear  combination  of  the  N  given  layer  values  .  The  S  has 

been  given  the  superscript  (v)  because  the  boundary  condition  on  vhas  been 
incorporated  in  the  solutions 

J  L  h«  to  be  differentiated.  L&  has  to  be  expressed  in  the  layer 
values.  This  can  be  done  as  follows:  From  eq.  (3^) 

Sr 


(  36  ) 


Substituting  eq.  (  35  )  into  eq  <  36  ,  end  interchanging  the  J  signs,  gives- 


( St )  =  L  vf  (r  l  s  ) 

__  r-  /  h 


(  37  ) 


So  the  derivative  of  the  dependent  variable  can  again  be  expressed 
combination  of  its  values  in  the  N  layers.  As  before  all  boundary  conditions 

have  been  incorporated  in  the  coefficients. 

Integration  gives  no  difficulty  either. 


Equation  (  34  gives: 

^  *  ,  r  ^  ^ 

cL<r  ^  4-  Vy 

u  k  iL 

o 


/ 


7T 
Aw 


Substituting  eq  (35  )  into  etl  )  gives: 

=  2L  y  'K  . 

4L  j.  L  *>L  t- 


J 


(  38  ) 


(  39  ) 


where 


<4—  «s-“ 


(  39a) 


it  should  he  pointed  out  that  these  coefficients  are  calculated  once  and  need  not 
be  calculated  again  unless  the  model  is  changed  to  a  different  number  of  layers. 

The  method  described  above,  has  been  applied  to  remove  differentiation 
and  integration  in  the  (Tcoordinate  in  the  equations  describing  the  Minta- 
Arakawa  model.  Though  the  latter  is  specifically  a  two -layer  model  (2  velocities,, 
the  vector  spherical-harmonies  model  will  not  be  limited  to  a  linearly  or 
quadratically  varying  velocity. 
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For  instance  the  hydrostatic  equation 


Now 


and 


X  ThLk)  Th 

k  1  *  ^<2 


is  the  boundary  condition  where  T^q  can  be  calculated  from  the  given 
ground  temperature 


Ts 

13  Lh 


Th 

lg 


As  before: 


■7“  to  (£) 

L  =* 


The  S(T)  are  calculated  by  inversion  of  the  matrix 
ki 

with  boundary  condition  eq  (  41  )  • 


of  equation  (  AO 


Substituting  eq.  (  AO  )  into  eq  (16a  )  results  in 


(16a  ) 

(AO  ) 

(41  ) 

(  A2  ) 

(  A3  ) 

(  44  ) 

(  45  ) 


with 
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(45a  ) 


There  are  N  independent  coefficients 
necessary. 


M(k) 


so  there  are  N  layers  values 


Suppose  that  originally  K  values  of  TM  are  given  as  the  initial  data.  K>  N. 
Now  the  ground  temperature  is  given  for  all  time.  So  equation  (  40  )  gives 
rise  to  K  equations  with  the  N  unknowns  T^(k)  (one  unknown  is  eliminated  with 
the  help  of  eq.  (  41  ))  Rather  than  integrating  K  layer  equations,  the  K  equations 


Th  T 
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(40a  ) 


are  solved  with  a  least  square  method,  giving 


r  - 1 

U  'ft  4  4-  l 


(46  ) 


N  new  "smoothed"  values  of  the  temperature  are  then  used  as  initial  values: 


4/  i-c.  kl  /  ' 


(47  ) 


TM  will  take  the  place  of  in  the  previous  discussion.  Also  the  T's  at 

L  j  L 

which  the  initial  data  are  given  need  not  be  the  same  as  then's  at  which  the 


integration  is  performed. 


Eq.  (  45  )  with  eq  (454)  takes  the  place  of  the  hydrostatic  equation. 
Next  the  equation  (17a  )  for  the  ground  pressure  will  be  transformed. 
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( 17a  ) 
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Again  A  is  written  as: 

'u 


(48  ) 


As  with  temperature,  V  is  subjected  to  boundary  conditions.  The  coefficients 

XM  must  be  calculated  by  inverting  eq  (  48),  adding  the  boundary  conditions 

L 

as  extra  equations.  An  example  of  this  will  be  given  later.  For  now  suppose 
the  solution  . is 


^  \u 


(49  ) 


where  the  index  i  runs  over  the  independent  layers 


Integrating  eq  (  48  )  gives 

i,  #(*)  j 
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Substitution  of  eq.  (  48  )  into  eq  (  50  )  results  in: 

/ 

AV-xxT  / 
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where 


(51a  ) 


The  equation  (  17a  )  for  the  ground  pressure  becomes: 


(52  ) 

<52a  ) 


Equation  (  49  )  resulted  from  adding  the  boundary  conditions  to  eq.  (  48)  and 
soiving  the  set  of  resulting  equations.  This  is  possible  because  eq  (  48  )  is  a 
system  of  N-l  equations  with  Nil  unknowns.  The  boundary  conditions^  two 
more  equation.  For  instance  in  the  case  of  a  given  drag  coefficient  d  =  £  ‘he 


boundary  condition  on  the  surface  is: 

,  /h7  \ 
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At  the  earth's  surface 


(/*=/,  so^is  not  a  function  of  (T\  Therefore 


1(H) 
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or  in  spherical  harmonics  coefficients: 


(  54  ) 


With  the  help  of  eq  (  ¥2  )  this  can  be  written  as : 


1 1  <  K7-  -L  % 


(  54^  ) 
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The  extra  equation  is: 


(  55  ') 


M(l). 


Together  .with  W-  0  at  (7*=  0,  which  gives  ^  ^ 


7>r 


This  forms  a  system  of  (N+l)' equations  ip  (N+l)  unknowns.  The  matrix 

i 


5  ^  1  is  the  inverse  of  the  coefficient  matrix  of  this  system, 

ki  J  <  •_  — 

In  this  case  there  are  (N-l)  layers  where  7  or  X  is  .calculated. 


is  represented 


by  a  cubic  equation 
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Let  £  =  — Q  =  .025.  The  four  equations  to  determine  the  XCL 


are: 
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f  M 


where  )(M  and  Xj  , 

4U 


are  measured  values. 


(  56  ) 


Inversion  of  the  matrix  in  this  equation  gives  the  coefficients  S ^  and  S^3 


where  k  =  0,  1,2,3.  (Also,  of  course  those  of  and  S^).  Dunpg  the 


;re 


computation  the  values  of  then's  on  the  ground  and  top  of  the  atmosphei 
need  not  be  used.  If  it  is  necessary  to  obtain  those  values,  they  can  be  calculated 

from  eq.  (  54  )  and  (54a  ) 
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where  theX/*^ 


are  found  from  eq.  (  55  ) 


Furthermore: 


Initially  four  or  more  7^.  might  have  been  given 

This  gives  rise  to  the  least  square  problem 


1  >  '<■  >  v 


to  be  solved  under  the  conditions  (54  )  and  (  55  ) 

If  the  solution  is  X  $  A  then  the  smoothed  values 

SL  ki 

to  be  used  in  eq.  (  56  )  are 

4 Li 


X 


Integration  of  eq.  (  H8a  )  goes  as  follows: 
From  eq.  (  51 ) 


where  they/^  are  given  by  eq.  (  51a  ) 


Further  more 
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It  follows  that  eq  (  18a  )  in  spherical  harmonics  coefficients  becomes: 
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where: 
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as 


Use  is  made  thatU?=0  for  (T -  0 

Also  it  should  be  kept  in  mind  that  if  is  still  a  function  of  (  However 

before,  if  TT  is  given  in  a  sufficient  number  of  points,  the  spherical  harmonics 

coefficients  00 


M  can  be  obtained  from  eq.  (  62  )  without  finite  differencing. 
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To  obtain  the  equation  of  motion  in  terms  of  spherical  harmonics 

coefficients,  three  terms  remain  to  be  evaluated,  - »  - 

2<r  a<r 


and  the  friction  term.  First  let 


4r\z. 

<U  L  4i 


Then 
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s  J?  and^ have  the  same  kind  of  boundary  conditions 
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Substituting  eq.  (  65  )  into  eq  (  64  )  gives: 
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(  66  ) 


So  the  derivitivea  are  again  a  linear  combination  of  the  (N-l)  independently 
calculated  layer  values.  In  the  same  way 
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where 


and 
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(  70  ) 
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Using  eqs(  66  ),  (  67  )•  (  68  ).  and  (  69  )  with  eqs  (  15a  )  and  (  15c  ) 


gives 


! 


(  71  ) 


In  a  similar  way  the  thermodynamic  equation  can  be  treated.  This  gives 
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With  oh  ,  *  Z>  k  r.S,.  where  is  calculated,  using 

Ji  k  ^  *'L 

the  proper  boundary  condition  on 

In  these  equations  terms  such  as  (  KW)  and  other  non-linear  terms 
are  to  be  calculated  with  the  transform  method. 

The  equation  contain  only  variables  in  the  (N-l)  layers, 

Eqs  (  71),  (  72  ),  (  73  )  together  with  eqs  (  52  )  have  to  be  solved.  They  are 

ordinary  coupled  first  order  non-linear  differential  equation.  The  boundary 
conditions  imposed  on  the  variables  &  and  ?  are  already  incorporated  in 
these  equations  and  they  have  to  be  solved  with  the  proper,  given,  initial 
conditions . 
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Angular  Momentum 


The  angular  momentum  vector  of  a  moving  fluid  of  density 
and  velocity  V  is  given  by 

dTT  (1) 

where  d'T'is  the  volume  element  and  r  is  the  radius  vector. 

To  begin  with  consider  the  case  of  a  sphere  surrounded  by  a 
fictitions  atmosphere  of  spherically  symmetric  density  • 

— 

V  can  be  expanded  in  the  form 


Remembering  that 
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The  contributions  for  these  three  kind  of  terms  are  then 


Consider  the  integral 
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M, 

Using  the  orthonormality  of  the  Y  *,,s  the  only  non -zero  terms  for  this 

Ij 

integral  must  be  those  originating  in  C^1  's  which  contain  terms.  A 

Xj  0 

look  at  equation  (  g  )  shows  that  there  are  no  such  terms. 
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Consequently  LD  =  0. 
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Again  only  terms  in  containing  will  contribute.  Equations  (6) 

^  Xj  0 

and  (8)  in  Section  1  show  that  these  terms  are 
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Thus  we  have  reached  the  result  that  to  the  approximation  that  p  zz  p  ) 
the  only  terms  capable  of  having  a  total  angular  momentum  different  from 
zero  are: 
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and 


) 
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Now,  of  course,  two  important  approximations  have  been  made:  the 
earth  was  assumed  spherical  and  the  density  was  assumed  to  depend  only 
on  the  distance  from  the  center  of  the  earth  but  not  on  the  angles  &  and  ^ > 
These  are,  however,  pretty  good  approximations  and  to  a  considerable 
degree,  then,  the  angular  momentum  of  the  atmosphere  resides  in  the 
three  terms  found  above.  In  particular  the  term  c!j*  (r,t)  C^(  ) 

is  the  main  component  of  the  zonal  circulation. 
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The  imaginary  value  does  not  matter,  it  simply  means  that  c^(r,t)  also 
has  to  be  imaginary  to  make  the  velocity  and  the  angular  momentum  real. 
The  result  implies  that  must  be  very  stable  since  if  we  assume  that 
at  one  instant  of  time  the  angular  momentum  only  has  a  z-  component  ^ 
i.  e.  c  j  (r,t)  =  c"j  (r,t)  =  0,and  the  other  velocity  terms  cannot  absorb 
any  total  angular  momentum,  the  law  of  conservation  angular  momentum 
will  keep  going.  The  next  step  is  to  remove  the  approximations 
and  to  write 
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Now 
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has  terms  of  the  form 
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Again  we  get  since  r  x  e^  =  0  that  =  0 


The  terms  L  i.  e. 
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We  now  see  that  unless  L-L'  we  get  no  contribution  to  L  also  unless 

B 

i 

M=M-1,  M'  or  M'+l  we  get  no  contribution.  Consequently  the 

calculation  becomes  relatively  simple.  For  L  one  gets  the  same 

c 

kind  of  result.  This  concludes  the  discussion  on  angular  momentum. 
Further  work  is  in  progress. 
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1.0  Introduction 


IBM  has  developed  analytic  methods  for  applying  Vector  Spherical 
Harmonics  (VSH)  to  the  global  weather  prediction  problem.  This  technique 
shows  promises  of  advantages  in  computational  speed  and  accuracy  over 
current  methods. 

The  Weather  Systems  Programming  department  is  currently  engaged  in 
developing  a  program  system  to  investigate  and  demonstrate  these  advantages. 
This  work  is  supported  by  IBM  internal  funding  (IRAD)  and  a  contract  with 
the  Advanced  Research  Projects  Agency  (ARPA) . 

The  IRAD  funds  cover  general  investigation  into  the  advantages  of  VSH 
methods,  while  the  ARPA  contract  requires  implementation  in  VSH  of  a  specific 
weather  model  -  the  Mintz-Arakawa  model.  A  general  program  system  has  been 
designed  to  address  both  of  these  objectives.  The  ARPA  contract  is  being 
charged  for  development  of  those  portions  of  the  system  necessary  to 
implement  the  Mintz-Arakawa  model.  IRAD  funds  are  being  used  for 
development  of  the  general  system.  The  ARPA  Contract  could  be  performed 
even  if  the  IRAD  task  were  terminated.  However,  the  performance  of  the 
ARPA  contract  will  be  significantly  enhanced  as  a  result  of  the  IRAD  task. 

It  may  be  of  interest  to  review  the  evolution  of  the  program  system. 

The  major  requirement  is  to  numerically  solve  models  expressed  in  VSH. 

This  involved  a  "predictor"  that  integrates  differential  equations  over 
time.  Accordingly,  the  contents  of  the  predictor  were  investigated  in 
detail.  These  included: 

o  Transformations  between  physical  and  spectral  space  including 
evaluation  of  differential  operators  in  physical  space. 

o  Application  of  numerical  integration  techniques  to  solve  the 
differential  equations  in  spectral  space. 

o  Preparation  of  an  initial  state  from  raw  data 

o  Output  capabilities 

u  Calculation  of  vector  coupled  sums  via  analytic  forms 

VSH  provides  a  method  for  solving  an  arbitrary  set  of  differential 
equations  in  spherical  coordinates.  It  was  desired  to  retain  this 
generality  in  the  predictor. 

Accordingly,  an  extremely  general  and  flexible  program  for  solving  an 
arbitrary  set  of  partial  differential  equations  in  spherical  coordinates 
(e.g.  a  weather  model)  was  designed  to  the  level  of  detail  needed  for  coding. 
This  program  is  highly  parameterized.  Choices  of  code  modules  to  be 
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included,  sizes  of  working  storage  arrays,  computational  loop  limits,  ^ 
data  set  format  specifications  in  the  job  control  language  can  all  b 
specified  via  parametric  expressions. 


This  report  discusses  the  following  topics: 

o  Project  requirements  that  influence  the  program  system  design 

o  Program  system  functional  requirements  &  design  objectives. 

o  A  development  plan  for  implementing  and  modifying  the  program 
system. 

o  The  current  state  of  the  program  system  design 


2.0  PROJECT  REQUIREMENTS  TO  BE  SATISFIED  USING 

PROGRAM  SYSTEMS 

2.1  Contractual  Requirements 

The  current  ARPA  contract  requires  a  test  computation  using  Vector 

Spherical  Harmonics  (VSH)  of  a  simplified  version  of  the  Mintz-Arakawa 
Model . 


2.2  Research  Requirements 

The  program  system  will  be  used  to  support  the  continuing  research  into 
VSH  methods,  properties  and  advantages  -  both  theoretical  and  compute-  ' 
tional.  The  following  sections  discuss  particular  research  areas  where 
programming  will  be  used,, 

2.2.1  Develop  Models 

VSH  permits  inclusion  of  expressions  in  the  differential  equations  ' 
(models)  hitherto  neglected  or  simplified  because  of  computational  com¬ 
plexity.  Solutions  of  models  incorporating  more  exact  and  complete 
descriptions  of  physical  phenomena  are  thus  possible.  More  accurate 
models  for  weather  prediction  and  climatology  can  be  developed  without 
being  unduly  constrained  by  computational  difficulties. 

2.2.2  Develop  Initialization  Techniques 

Weather  data  collected  on  a  global  scale  is  not  of  consistent  accuracy  and 
completeness.  This  is  particularly  true  of  wind  data.  This  inconsistency 
is  further  compounded  by  the  nonsynchronous  nature  of  satellite  collected 
data.  Initialization  of  a  global  predictor  requires  resolution  of  this 
data  into  a  complete  sample  of  greatest  possible  accuracy.  VSH  techniques 
will  be  developed  to  produce  this  initial  condition. 


2.2.3  Eliminate  Radial  Dependence 

Application  of  VSH  to  partial  differential  equations  (PDE ' a)  in  spherical 
coordinates  produces  PDE's  with  radius  and  time  as  the  independent 

variables  and  VSH  expansion  coefficients  118  5h® J^i^depen- 
Methodology  will  be  developed  to  transform  these PDE  s  ini t  ti  A  ^ 

dent  ordinary  differential  equations.  The 

by  expanding  each  VSH  coefficient  in  suitable  radial  functions. 


2.2.4  Develop  Computational  Efficiency 

Application  of  VSH  to  a  complex  predictive  model  produces  a  tractable 
«d  computationally  efficient  formulation.  Algorithms  employed  in 
su  h  as  transformations  and  integration  as  well  as  da“  ““**“*“ 

SqueO  will  be  developed  to  further  enhance  the  computational  efficiency 

of  the  predictions. 


2.2.5  Determine  Accuracy 

Prediction  accuracy  end  conpot.tlon.1  speed  depend  on  .  conblnatlon  of: 

a.  The  model  used 

The  quality  of  the  initial  conditions 
The  truncation  limits  of  VSH  expansions 
The  form  and  order  of  radial  expansions 
e.  The  computational  algorithms  employed 

The  sensitivity  of  predictive  ^\?Zie5£d  “ 

well  as  the  effects  of  coupling  among  these  areas,  will 


b. 

-s 

C. 


d. 


3.0  PROGRAM  SYSTEM  FUNCTIONAL  REQUIREMENTS 


3. 1  Integration 

The  major  function  of  the  program  system  Is  to  numerically  solve  the 
model  equations  using  VSH  methods.  This  Involves  a  prediction  that 
generally 

a.  Expands  physical  variables  In  vector  and  scalar  spheri¬ 
cal  harmonic  coefficients  and  transforms  harmonic 
coefficients  to  physical  variables 

b.  Computes  linear  and  nonlinear  terms  in  differential 
equations  to  form  derivatives 

c.  Integrates  over  time  in  VSH  coefficients  to  produce 
predictive  results. 

The  exact  method  of  integration  to  use  is  not  known  at  this  time.  One  of 
the  research  objectives  is  to  find  the  best  applicable  integration 
method.  Numerical  integration  methods  are  generally  divided  into  two 
classes: 

a.  Explicit,  or  open,  methods  in  which  the  state  at  a  given 
time  is  computed  as  a  linear  combination  of  state  and 
derivative  values  at  previous  time  points. 

b.  implicit,  or  closed,  methods  which  are  the  same  as 
explicit  methods  with  the  addition  of  some  multiple 

of  an  estimate  of  the  derivative  evaluated  at  the  same 
time  point  at  which  the  state  is  being  computed. 

The  integration  method  may  be  a  single  step,  where  only  the  state  and 
derivative  values  from  one  past  time  are  needed,  or  a  multi-step,  where 
state  and  derivative  values  from  several  previous  times  are  required. 

A  requirement  for  integration  is  the  computation  of  derivatives  in  spec¬ 
tral  space. 

The  elements  of  spectral  derivative  computation  are: 

a.  Computation  of  radial  derivatives  of  spectral  variables, 
which  is  itself  a  research  area. 

b.  Formation  of  differential  terms  in  physical  space  from 
spectral  variables  and  differential  terms. 

c.  Transformation  of  the  differential  equations  in  physi¬ 
cal  space  to  corresponding  differential  equations  in 
spectral  space* 


3.2  Initialization 

The  predictor  requires  an  initial  state  of  ^  ^ 

for  Its  nroenostic  computations.  Derivation  of  this  initial  stare 
r»  u  functional  requirement  of  tha  progra.mlng.yeem 

Current* technique,  are  generally  regarded  “ 

LT 

SSSETaS  mixes  of  raw  data.  Thla  nay  be  done  In  several  atages. 


a.  Artificially  generate  an  initial  state 

b.  Apply  VSH  to  existing  initialization  methods 

c.  Implement  currently  unknown  VSH  methods. 


3. 3  Output 

Data  fro.  the  predictor  must  be  output  to  Inform  ‘’“™'£^1L“He“’ 

for° the'predlc torque t'b^output^Furth.r,  benchmark  comparisons  and 
contractual  requirements  nay  necea.ltat.  particular  content  and  format 
of  output. 


i.4  Table  Generation 

lata  used  In  the  prediction.  Initialisation  and  output  will,  whenever 
•oeelble,  be  precomputed  to  enhance  computational  speed. 

b.  a  la„c  .mount  of  data  in  the  tables  (possibly  several 
ilssion  quantities)  and  computing  them  is  a  major  task.  The  ta  £« 

predictor. 


3.5  Flexibility 

system,  from  the  highest  level  to  the  lowest. 
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3.5.1  Model  Changes 


Models  are  sets  of  PDE's  in  spherical  coordinates  and  their  boundary 
conditions.  The  program  system  must  be  able  to  accept  the  following  types 
of  model  changes  which  are  enumerated  in  order  of  complexity. 

a.  Execution  time  parameter  changes  (e.g.,  assign  values  to 
the  symbolic  coefficients  that  appear  in  the  model 
equations) 

b.  Reformulation  of  an  existing  equation  (e.g.,  addition  of 
terms  to  more  closely  approximate  the  physical  system) 

c.  Modification  of  boundary  conditions  (e.g.,  change  a 
constant  boundary  condition  to  a  function  of  time) 

d.  Add  equations  (e.g.,  explicitly  model  the  dynamics  of  a 
physical  variable  hitherto  assumed  constant) 

e.  Change  physical  system  (e.g.,  model  seismo logical 
phenomena). 


3.5.2  Algorithm  Changes 

The  computational  algorithms  used  to  implement  the  various  transformations, 
differential  equation  component  calculations  and  integration  are  themselves 
the  subject  of  investigation  as  to  accuracy  and  efficiency. 

Possibly  several  candidate  algorithms  will  be  tried  in  a  given  application 
and  modifications  to  specific  algorithms  will  be  made.  This  implies  the 
need  for  ease  of  replacement  and  modification  of  algorithmic  components  of 
the  programming  system. 


3.5.3  Experimental  Changes 

The  VSH  method  expands  physical  variables  in  spectral  functions  of  latitude 
and  longitude  and  currently  unknown  functions  of  radius.  The  program  sys¬ 
tem  must  support  experimental  changes  of  the  following  forms  in  these 
expansions : 

a.  Number  of  terms  in  spectral  expansion 

b.  Functions  to  be  used  in  radial  expansion 

c.  Number  of  terms  in  radial  expansion. 
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3.5.4  Initialization  Changes 

There  are  three  main  initialization  objectives 

a.  Apply  VSH  to  efficiently  produce  an  initial  sample  of 
greatest  possible  consistency  and  accuracy.  This  sample 
will  contain  the  variables  necessary  to  initialize  our 
model. 

b.  Accept  a  verity  of  in^ut  data  and  formats  and  transform 
them  to  the  iorm  necessary  to  initialize  our  model.  This 

is  particularly  true  when  performing  benchmark  comparisons. 

c.  Experiment  with  the  quantity  of  initialization  data  neces¬ 
sary  for  accurate  prediction. 

These  objectives  imply  continuing  changes  and  additions  in  the  Initializa¬ 
tion  section  of  the  program  system. 


3.5.5  Output  Formats 
There  are  two  types  of  output: 

a.  General  predictive  results  describing  the  weather  state. 
The  form  and  content  is  generally  well  defined. 

b.  Data  t:o  illuminate  particular  points  of  interest  in 
research.  The  form  and  content  is  usually  defined  - 
and  redefined  -  as  experimentation  proceeds. 

Also,  the  content  and  format  of  output  will  be  determined  by 

a.  Internal  Use 


b .  Contract  Requirements 

c.  Benchmark  Comparisons. 

Varieties  of  data  must  be  aggregated  and  output  in  desired  formats;  plots 
of  various  types  must  be  produced.  The  data  output  may  be  raw  or  processed 
input  data,  tables  of  coefficients  used  in  the  computations,  the  results 
of  intermediate  steps  in  the  predictions,  or  the  predictions.  In  many 
cases  the  exact  nature  of  the  output  will  be  defined  after  the  baseline 
system  structure  has  been  built  (changes  and  additions  will  be  required 
to  satisfy  previously  undefined  or  unanticipated  requirements).  The  out¬ 
put  section  must  be  able  to  satisfy  these  data  needs,  especially  those 
generated  by  specific  research  questions,  contract  requirements  and 
benchmarks . 
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4.0  PROGRAM  SYSTEM  OBJECTIVES 


The  program  system  must  be  developed  to  meet  two  major  objectives 

a.  To  provide  a  general  and  flexible  testbed  for  investi¬ 
gation  in  the  major  research  areas. 

b.  To  provide,  when  necessary,  the  most  efficient  overall 
program  for  the  solution  of  a  given  problem. 

4.1  Baseline  and  Particular  Systems 


4.1.1  Baseline  System 

The  first  objective  leads  to  a  more  elaborate  and  perhaps  less  efficient 
overall  structure  than  might  be  attainable  if  certain  segments  of  the 
system  could  be  made  less  general.  However,  a  general  structure  is  nec¬ 
essary  since  the  major  areas  for  research  impact  every  aspect  of  the 
program  system  and  require  maximum  flexibility  for  satisfactory  support 
and  response.  The  format  of  the  model  equations,  the  boundary  conditions, 
the  radial  fitting  functions  and  the  initializing  process  are  all  certain 
to  be  changed  several  times  in  the  search  for  maximum  accuracy  and 
fidelity  to  physical  reality.  Algorithms  and  program  parameters  for  such 
things  as  transformations,  derivative  calculations  and  integration,  are 
bound  to  be  altered  and  replaced  to  enhance  accuracy  and  computational 
efficiency.  Therefore  a  system  must  be  designed  to  accept  changes,  both 
trivial  (integration  time  step)  and  fundamental  (integration  method) 
without  undue  effort  and  delay.  We  will  call  this  a  "Baseline  System,' 
in  that  it  provides  the  framework  in  which  various  alternatives  can  be 
easily  implemented  and  investigated. 

The  Baseline  System  must  provide  a  flexible  and  easily  changeable  testbed 
for  experimentation  in  the  major  research  areas.  It  must  serve  as  the 
development  ground  and  checkpoint  for  particular  systems. 


4.1.2  Particular  Systems 

Generality  serves  the  overall  research  objectives,  but  there  may  be 
occasions — baseline  comparisons,  contract  requirements,  promising 
research  configurations  when  the  maximum  overall  program  efficiency  is 
desirable. 

The  information  gained  from  research  using  the  baseline  system  will  allow 
a  given  configuration  of  model  equations  and  computational  algorithms 
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4.3  Computer  Requirements 

StLSSSCJISS:  tTT  en0r"°“8  a”°“"“  of  *».  »-th  initial 
for  eralSr  Furthf »  data  requirements  are  dynamic  depending, 
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such  a  system  large  machines  with  sufficient  core  capacity  for 

cult.  A  better  obi erMt/  access  to  such  a  machine  may  be  diffi- 
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4.4  Parameterization 

The  program  system  development  must  account  for  implication  of  the  nrn 

o«p“:1s.^^nr^:;.values•  i"te8ration  ata»  — u*t,  .„d 

“"cuMon^.rH  lmPJlM  that  tL'  aB0Unt  0f  “ra  squired  for  a  given 
xecution  will  depend  on  execution  time  parameters  such  as  expansion 

SSTSuST;  alaaa  °f  -««•  rnb'^s  5n  va“So- 

thms  will  also  depenc  on  these  execution  time  parameters.  8 

Flexibility  in  Models  and  algorithms  and  generality  of  innut  a„a  «  * 
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In  summary,  there  are  a  large  number  of  parameters  necessary  to  define 
a  given  system  configuration.  They  are: 

a*  Execution  time  parameters,  such  as  spectral  expansion  trun- 
cation  limits ,  some  of  which  determine  compile  time 
parameters 

b.  Compile  time  parameters,  such  e.s  program  segments  and  core 
segmentation,  some  of  which  are  functions  of  execution 
time  parameters. 

c.  Data  sr.t  definition  parameters,  such  as  record  length, 
some  of  which  are  functions  of  execution  time  parameters. 

It  is  desirable  to  describe  a  given  system  configuration  by  specifying  as 
few  of  these  parameters  as  possible.  This  can  be  done  by  determining 
those  parameters  that  are  independent  and  used  in  specifying  other 
parameters.  It  should  be  necessary  to  specify  only  those  parameters. 
Algorithms  should  be  built  into  the  system  to  determine  values  for  the 
other  dependent  parameters,  and  to  construct  the  code  segments  and  JCL 
automatically . 

This  cau  be  done  through  the  use  of  precompilation  facilities  described 
in  the  Appendix.  Briefly,  this  involves  the  PL/1  precompilation  facility 
operating  on  FORTRAN  and  JCL  —  as  well  as  PL/1— code.  Parametric  expres¬ 
sions  in  the  code  are  evaluated  to  determine: 

a.  Segments  of  code  to  be  included  in  the  compilation. 

b.  Storage  allocations. 

c.  DCB  specifications  in  the  JCL. 

This  implies  that  the  system  must  be  recompiled  when  execution  time  param¬ 
eters  which  affect  precompilation  parameters  change. 

Also,  the  iterative  nature  of  the  research  and  the  number  of  parameters  in 
the  system  add  an  additional  constraint.  A  particular  version  of  the  sys¬ 
tem,  with  its  execution  and  precompilation  parameters  set,  must  be 
archived  and  capable  of  easy  restoration. 
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5.0  BASELINE  SYSTEM  DEVELOPMENT  PLAN 


Because  of  its  gene^^lity  and  s izz,  a  coherent  plan  for  the  development 
of  the  Baseline  System  is  necessary.  The  plan  and  the  implementation 
procedure  should  allow  the  system  structure  to  be  visible  and  easily 
modifiable  during  its  construction.  The  system  should  be  built  so  that 
changes  at  one  stage  of  the  implementation  do  not  cause  extensive  modi¬ 
fication  of  previously  written  code. 


5.1  Structured  Programming 

Baseline  System  design  and  development  will  use  a  technique  called 
structured  programming  in  order  to  maintain  visibility,  consistency  and 
efficiency — both  in  implementation  and  execution.  The  basic  tenet  of 
structured  programming  is  that  computer  programs  can  be  designed  and 
written  with  a  high  degree  of  structure,  which  permits  them  to  be  more 
easily  understood  for  testing,  maintenance,  and  modification. 

To  further  enhance  the  readability,  visibility  and  coherence  of  the  code, 
a  repetitive  process  of  code  segmentation  is  undertaken.  The  first  step 
in  this  process  is  to  formulate  a  one-page  skeleton  program  which  repre¬ 
sents  the  entire  program.  This  is  done  by  selecting  some  of  the  most 
important  lines  of  code  in  the  program  and  then  filling  in  what  lies 
between  those  lines  by  names.  Each  new  name  will  refer  to  a  new  segment 
to  be  stored  in  a  library  and  called  by  a  macro  facility.  In  this  way, 
we  produce  a  program  segment  with  something  under  50  lines,  so  that  it 
will  fit  on  one  page . 

Now,  the  segments  at  the  next  level  can  be  written  in  the  same  way,  refer¬ 
ring  as  appropriate  to  segments  to  be  later  written  (also  setting  up 
dummy  segments  as  they  are  named  in  the  library).  As  each  dummy  segment 
becomes  filled  in  with  its  code  in  the  library,  the  recompilation  of  the 
segment  that  includes  it  will  automatically  produce  new  updated,  expanded 
versions  of  the  developing  program. 

This  leads  to  a  top  down  system  implementation.  The  system  is  created  in 
execution  sequence.  Only  instructions  which  can  be  executed  as  they  are 
created  are  coded,  all  prior  instructions  required  have  already  been 
coded.  In  the  top-down  approach,  the  highest  level  code  is  written  first:, 
then  the  next  lower  level  code,  etc.,  down  to  the  lowest  levels  of  code. 

Under  OS/360,  for  example  job  control,  linkage  editor,  "supervisory"  and 
"data  management"  code  written  in  that  order,  and  only  then  the 
source  modules  which  typically  give  a  system  its  functional  capability. 
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Thus,  development  proceeds  through  the  controlled  addition  of  new  modules 
to  an  always  checked-ouc  program.  That  Is,  supervisory  programs  run  early 
in  the  development  phase,  first  calling  on  dummy  modules  for  which  func¬ 
tional  modules  are  substituted  later. 

To  augment  the  ordered  implementation  of  the  top-down  approach,  record 
keeping  and  program  module  library  procedures  .mist  be  formatted  and 
followed.  In  the  structured  programming  approach,  this  is  done  via  the 
Programming  Production  Library  (PPL). 

The  PPL  is  used  to  maintain  the  current  status  and  past  history  of  all 
code  generated  in  a  project  in  a  human  readable  form,  except  as  specifi¬ 
cally  noted  to  the  contrary.  The  current  status  of  the  PPL  will  include, 
at  least,  separate  Library  Sections  which  represent  the  following  six 
classes  of  code  and  data: 

a.  Job  Control  Source  Code 

b.  Linkage  Editor  Source  Code 

c.  Program  Source  Code 

d.  Object  Module  Code 

e.  Load  Module  Code 

f.  Test  and  Control  Data 

These  libraries  generally  reside  in  partitioned  data  sets  on  direct  access 
storage  devices. 

Standard  PPL  procedures  exist  to  perform  the  following  functions: 

a.  Add/Delete/Modify  any  Source  Member 

b.  Assemble/Compile  any  Object  Member 

c.  Linked! t  any  Load  Member 

d.  Execute  Debug/Test  Program. 


5.2  Subsystem  Structure 

In  order  to  best  satisfy  the  research  requirements  and  program  system 
functional  requirements  and  objectives,  the  baseline  system  will  be  seg¬ 
mented  into  three  subsystems. 


5.2.1  The  SO  Subsystem 
This  subsystem  will  include 

a.  Integration  of  the  model  (prediction) 

b.  Initialization  for  prediction 

c.  Output 

d.  Table  generation  for  prediction,  initialization,  and 
output . 

Each  of  the  above  functions  -  Initialization,  tabic  generation,  predic¬ 
tion,  output  -  addresses  various  project  requires  ants .  Each  function 
contains  many  options  and  the  functions  may  be  combined  in  various  ways 
depending  on  what  is  desired.  The  GO  subsystem  will  therefore  consist 
of  each  function  implemented  as  a  separate  program  module.  The  execution 
of  each  function  in  turn,  and  the  options  chosen  in  each  function,  will 
be  controlled  by  parameters  specified  by  the  user,  or  constructed  from 
user  specified  parameters. 


5.2.2  The  BUILD  Subsystem 

The  purpose  of  this  subsystem  is  to  construct  an  executing  version  of 
the  GO  subsystem.  It  automates  to  the  greatest  extent  possible  the  high 
degree  of  parameterization  present  in  the  GO  subsystem. 

The  BUILD  subsystem  will  do  the  following: 

a.  Use  precompilation  facilities  to  Include  code  segments  and 
arrange  working  storage. 

b.  Construct  JCL  for  executing  the  GO  subsystem  and  insert  it  in 
Procedure  Library. 

c.  Compile  and  Link  edit  Input,  Table  Generation,  Predictor, 
Output  load  modules. 

d.  Construct  control  parameter  data  records  to  select  options 
within  the  modules  and  choose  the  modules  to  be  executed. 

e.  Archive  the  GO  subsystem  including  source  and  object  code, 
load  modules,  JCL  and  control  parameter  data  records. 


01; 
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5.2. J  The  RESTORE  Subsystem 


The  purpose  of  this  subsystem  is  to  restore  a  GO 
complete  with  JCL  and  control  parameter  records, 
the  BUILD  subsystem.  Only  that  portion  of  the  GO 
execution  will  be  restored. 


subtfystfsm  configuration, 
that  was  archived  by 
subsystem  necessary  for 


The  following  will  be  done: 


a’  Restore  the  Input,  Table  Generator,  Predictor  and  Output  load 
modules  to  the  direct  access  load  library. 

b.  Restore  the  GO  JCL  procedure  to  the  procedure  library 
for  subsequent  execution. 


Restore  the  control  parameter  data  records  to 
parameter  direct  access  data  set. 


the  control 


5.3  Subsystem  Execution  Control 

Execution  control  within  a  subsystem  and  communication  between  the  sub¬ 
systems  are  implemented  by  the  Control  Parameter  Preparation  (CPP>  it 
gram  and  an  associated  direct  access  data  set  -  the  ££  5a™  sit.  ? 

The  CPP  data  set  consists  of  various  types  of  records: 


a. 


System  Master  Record  -  Major  record  for  subsystem 
munlcatlon,  contains: 


com- 


1.  A  code  identifying  the  system  configuration  for 
archiving  purposes. 


b. 


Indicators  for  GO  subsystem  module  execution  inclusion 

3.  Pointers  to  compile  time  parameter  records  to  be  moved 
into  the  subsystem  master  records  at  build  time 


^t6r  Records  •  °“e  each  for  Input,  Predict, 
and  Output  modules.  Identifies  all  requisite  parameters  for 
each  function. 


Compile  time  parameters 


Pointers  to  data  records  containing  execution  time 
parameters. 


2. 


c.  Parameter  Records  -  Contain  either  execution  or  compile 
time  parameters  for  the  GO  subsystem  modules. 

d.  Working  Records  -  Used  by  CPP  program  for  CPP  file 
maintenance . 

at  d  rha  CPP  nroeram  is  to  maintain  the  CPP  file.  Data 

functions  are  performed  in  the  BUILD  and  CO  subsystems. 

In  addition,  the  tasks  performed  by  the  CPP  program  in  each  of  the  sub- 
systems  are: 


b. 


GO 

1. 

2. 

3. 


Checks  CPP  master  and  data  records  for  completeness 
and  validity 

Sets  an  operating  system  condition  code  to  control 
execution  of  desired  GO  Subsystem  modules. 

Marks  CPP  master  record  to  indicate  that  tables 
have  been  generated  or  copied. 


BUILD 
1. 

2. 

3. 

A. 

5. 


Checks  input  parameters  for  validity  and  complete- 
ness 

Updates  the  CPP  Master  record  according  to  input 
parameters 

Marks  the  system  master  record  "No  Go"  until  BUILD 
is  complete 

Marks  the  CPP  master  record  for  table  generation 
in  the  GO  Subsystem 

Constructs  an  input  data  set  of  utility  control 
statements  for  the  archiving  step. 


restore 
1. 


2. 


Marks  the  system  master  record  as  "No  Go"  until 
restoration  is  complete 

Constructs  an  input  data  set  of  utility  control 
statements  for  the  step  that  restores  from  the 
archive  data  set 
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3.  Harks  the  system  master  record  for  the  GO  Sub¬ 
system  to  copy  tables  from  the  archive  data  set. 

4.  Updates  the  CPP  master  record  from  the  archive 
tape. 


5.4  Baseline  System  Development  Procedure 

Figure  5-1  presents  the  procedure  by  which  the  Baseline  System  will  be 
developed.  Mechanisms  are  provided  for  system  generation  and  subsequent 

revisions . 

There  is  a  three  level  classification  for  Baseline  Systems: 


a.  Version 

Recompilation  with  different  compile  time 
parameters.  No  new  source  code,  but  possible 
corrections  to  existing  code 

b .  Extension 

Modification  addition  of  new  options  (additive 
growth).  The  key  is  a  CPP  program  addition  to 
recognize  and  manage  the  new  options. 

c.  System 

A  new  system  has  been  generated  when  the  following 
cannot  be  done 

3 .  Allow  any  restored  version  to  execute 

2.  Rebuild  any  earlier  version  with  suitable 
input  parameters 
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Set  system  number,  extension  number  and  version 
number  to  Zero 


End  of  Project 


Figure  5-1.  System  Development  Procedure 
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6.0  BASELINE  SYSTEM  DESIGN 


This  section  describes  in  detail  the  three  subsystems  of  the  Baseline: 
System.  The  CPP  program  as  described  in  Section  5.4  is  referenced  in 
each  subsystem  discussion.  The  data  sets  and  their  use  ivi  each 
subsystem  is  summarized. 


6.1  The  GO  Subsystem 

The  major  function  of  the  GO  Subsystem  is  to  produce  predictive  estimates 
of  the  weather.  This  involves  four  functions: 

a.  Table  Generation  -  Compute  and  write  data  sets  of  constants 
used  in  the  GO  Subsystem.  These  data  sets  may  optionally  be 
copied  from  archived  data.  This  function  ic  included  in  the 
GO  Subsystem,  rather  than  the  BUILD  (or  RESTORE)  subsystem, 
because  the  JCL  describing  the  table  data  sets  depends  on 
precompilation  parameters.  This  JCL  is  built  as  part  of 
the  GO  procedure. 

b.  Input  -  Produce  an  initial  sample  for  prediction  from  raw 
data.  This  may  be  done  in  several  ways: 

1.  Artificially  generate  an  initial  state 

2.  Apply  VSH  to  existing  initialization  methods 

3.  Implement  currently  unknown  VSH  methods 

c.  Prediction  -  Integrate  the  model  equations.  This  involves 
numerical  Integration  of  the  model  equations  expressed 

in  VSH.  The  Integration  method  is  the  subject  of  research. 
Derivatives  in  VSH  must  also  be  computed  for  the  integra¬ 
tion. 

d.  Output  -  Produce  output  concerning  any  of  the  above  func¬ 
tions  .  Tailor  output  to  meet  benchmark  specifications  and 
contract  requirements. 

The  GO  Subsystem  will  consist  of  ioad  modules  for  each  of  the  above  func¬ 
tions  as  separate  job  steps.  These  load  modules  will  be  produced  by  the 
BUILD  Subsystem  discussed  in  Section  6.2. 

There  are  many  execution  options  possible  among  the  four  functions.  For 
example,  the  Input  module  might  process  some  raw  data  and  pass  the  results 
to  the  Output  module  for  plotting.  The  Table  Generator  and  Predictor  mod¬ 
ules  are  not  executed. 


The  execution  of  a  given  function  will  be  controlled  by  operating  system 
condition  codes .  The  CPP  program  will  set  the  condition  code  indicating 
execution  of  any  combination  of  the  modules .  The  modules  desired  will  be 
indicated  via  input.  CPP  will  also  perform  a  validity  and  completeness 
check  on  the  CPP  data  set  after  it  has  processed  any  input  parameters. 

If  validation  fails,  CPP  sets  a  job  condition  code  of  99,  which  terminates 
execution.  The  functional  flow  of  the  GO  Subsystem  with  attendant  data 
sets  is  shown  in  Figure  6-1.  The  data  sets  are  discussed  in  the  following 
section,  and  the  functions  are  discussed  in  detail  in  subsequent  sections. 


6.1.1  GO  Data  Sets 


6. 1.1.1  TABLES 

Precomputed  constants  are  used  in  the  Input,  Predictor  and  Output  modules. 

They  reside  on  a  number  of  TABLES  data  sets  constructed  by  the  Table 

Generation  Module.  Unless  otherwise  noted  they  are  direct  data  sets. 

a.  TABLE  CONTROL  -  Sequential  data  set  containing  copies 
of  the  CPP  Master  records  at  table  construction  time. 

b.  PREDINIT  -  Sequential  data  set  containing  various  con¬ 
stants  that  will  reside  in  core  during  Predictor  Module 
execution.  They  are  tables  defining  radial  fitting  func¬ 
tions,  frequently  used  Clebsh-Gordan  coefficients  and 
constants  used  in  the  Fourier  step  of  the  transformation 
from  spectral  to  physical  space. 

c.  RADIAL  ONE  -  Constants  used  in  constructing  first  order  radial 
differential  operators  of  spectral  coefficients. 

d.  RADIAL  TWO  -  Constants  used  in  constructing  second  order  radial 
differential  operators  of  spectral  coefficients. 

e.  PS  LEGENDRE  -  Constants  used  in  the  Legendre  step  of  the 
transformation  from  physical  to  spectral  space. 

f.  PS  FOURIER  -  Constants  used  in  the  Fourier  step  of  the 
transformation  from  physical  to  spectral  space. 

g.  VS  REAL  -  Constants  used  in  transforming  between  VSH  and  Scalar 
Spherical  Harmonic  (SSH)  real  components. 

h.  VS  COMPLEX  -  Constants  used  in  transforming  between  VSH 
and  SSH  complex  components 

i.  SP  LEGENDRE  -  Constants  used  in  the  Legendre  step  of  the 
transformation  from  Spectral  to  physical  space. 
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6. 1.1. 2  ARCHIVE 


TablenGe^rftta  ^°ntains  a  C0Py  of  the  TABLE  data  sets.  The 

Table  Generation  Module  writes  this  file  when  it  builds  tables,  and 

reconstructs  them  from  this  file  when  a  system  has  been  llstlled 


6. 1.1. 3  RAW  DATA 

dVa  aet0  containing  weather  obaervatlons .  The  for.  and  con- 

ziu. not  currentiy  known-  *“•  »iu  be  «*•  »«*-* 


6. 1.1. 4  IM  OUT 

set*  containing  the  initial  conditions  needed  by  the 
Predictor  Module.  These  conditions  will  be  expressed  in  physical  space 

pr»rdVr::;.Ihe  forMt  °f  this  -« *  «■«  «  -tk 


6.1. 1.5  PM  IN 


Sequentisi  data  set,  containing  initial  conditions 
Module.  This  is  either  a  previously  computed  IM  OUT 


needed  by  the  Predictor 
or  PM  OUT  data  set. 


6. 1.1. 6  PM  OUT 

Sequential  data  set  containing  Predictor  Module  output.  It  consists  of 

The  foraatVofUthi°fHthe  phW8JCal  and  8Pectral  variables  at  selected  times 
The  format  of  this  data  set  is  the  same  as  that  of  the  IM  OUT  data  set. 


6. 1.1. 7  OM  IN 

!™:Ldata  8et’  4containin8  data  to  be  processed  by  the  Output  Module. 
This  is  either  a  previously  computed  IM  OUT  or  PM  OUT  data  set. 


6.1. 1.8  OM  OUT 

?°re  data  8ets  c°ntaining  GO  Subsystem  results  in  user  specified 
S? table8»  pl°t8'  tapes,  etc.).  The  form  and  content  of  these 
W  Subs^t^.n  CUrrently  known*  11118  wil1  be  the  primary  output  of  the 
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,  I  ■ 

6. 1.1, 9  Work  Data  Sets 

argali  direct  P'edic£or  M°dule  storage  without  undue  use  of  core.  They 
execution  data  sets  and  exist  only  during  Predictor 

an^its^pectra^coeff icients^They  ^  Phy8±Cal  Pf°~  ™iabla 

i  i  i  ■  *  1 

a.  P  VALUE  -  The  value  of  the  variab.’e  in  physical  space. 

b.  P  DIFF  -  The  va^ue  of  spatial  differential  operators 

appxxed  to  the  physical  variable.  i 

c.  S  VALUE  -  The  value  of  the  variable  in  spectral  space  at 
one  or  more  time  instants. 

,  1  | 

d.  SDR  -  The  value  of  the  first  and  second  radial  partial 
derivatives  of  the  spectral  coefficients  for  this 
variable.  1 

i  i 

i 

e.  SDT  -  The  value  of  the  temporal  Partial  derivatives  of 

the  spectral  coefficients  for  this  variable  at  one  or 
more  time  instants.  1 

1  1  '  i 

SrLd“l?l8°  0“  ““  SCt  f°r  311  the  variables  Indicated  In 

:  1  il 

spacf?  ~  TC  VSlUe  °f  thC  dld8no8tic  variables  in  physical 

\  * 

6.1.2  GO  Subsystem  Components  i  , 

,  1  '  :  . 

6.1.2. 1  Table  Generation  Module 

i  I  i 

This  module  performs  one  of  two  functions  depending  on  the  status  of  th» 
table  flag  iu  -he  System  Master  Record  of  the  CPP  data  set: 

fl‘  SS**  W^,te  th£  TABLE  data  8et8  lf  ^ble  build  is 

the  Sramfdat!  s"!!"”1"68  arChl'rInS  tha  IableB  °r 

b‘  tab le'res to ration  **  "e  “  , 
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A  detailed  design  of  thl.3  module  has  not  yet  been  done.  However,  the 
format  and  contents  of  the  tables  are  known,  as  are  the  formulae  neces¬ 
sary  to  produce  them.  Implementation  of  this  module  will  be  in  the  most 
straightforward  fashion,  since  use  and  not  production  of  the  tables  is  of 
prime  interest. 


6 . 1 . 2 . 2  Input  Module 

This  module  produces  an  initial  sample  for  the  predictor  from  raw  data 
using  VSH  methods.  Input  to  the  module  is  the  IM  IN  data  set  and  module 
output  is  the  IN  OUT  data  set.  Of  the  three  major  functions  of  the  input 
module  -  input,  processing,  output  -  only  the  format  and  contents  of  the 
output  data  set  are  currently  known  in  any  detail. 

The  raw  data  may  come  in  several  forms  and  there  may  be  varying  mixes  of 
data.  This  is  particularly  true  in  the  case  of  benchmarks  and  implementa¬ 
tion  contracts,  where  the  form  and  content  of  the  raw  data  will  be 
specified  as  a  part  of  the  contract. 

Implementation  of  input  processing  is  currently  seen  as  a  three  stage 
process. 

a.  Artificially  generate  an  initial  state  to  begin  exercising 
VSH  predictive  formulations. 

b.  Apply  VSH  methods  to  existing  initialization  technology. 

c.  Research  and  implement  currently  unknown  VSH  methods. 

As  in  Table  Generation,  implementation  will  be  straightforward.  The  princi¬ 
pal  merit  of  VSH  methods  in  this  area  is  a  gain  in  the  consistency  and 
accuracy  of  the  initial  sample.  The  proportion  of  time  spent  in  initiali¬ 
zation  relative  to  prediction  obviates  undue  concern  for  coding  efficiency 
here . 


6. 1.2. 3  Predictor  Module 

The  prime  function  of  the  Predictor  module  is  to  integrate  the  model 
equations  expressed  in  VSH  coefficients.  This  involves  an  initialization 
consisting  of 

a.  Reading  and  setting  up  incore  tables 

b.  Building  first  copies  of  applicable  work  data  sets  from 
initial  data.  This  involves  data  sets  such  as  P  Value 
.ad  S  Value . 


c. 


Writing  the  initial  state  as  the 
PM  OUT  data  set. 


first  record  of  the 


I , "  1  “ ra “vc  integration  process  is  entered.  Startle* 


^f|“r  ”0d^e  Snd  its  attendant  data  sets  have  b»en  designed  with 

S*"i\h^^^ 

dfriwtivf 1"VOlVe  ,0rmln8  llnear  “"bination,  of  past  states 

r'ftSSaSasSSSsS?1 

rSSSFdv-?—  ' Sri  £  date 

^ectra'/sp’ale'^TMs^.  18  the  “"P“tntion  of  derivatives  in 

formation  o?  the^ne^r  e™M  ??  "f"  °f  the  integration  code.  The 

auction  or  tne  linear  combination  of  Btates  and  derlvat-iwa  f„, 

.  “d  Cha"8ad  Separately  from  the  derivative  coilpota- 

actuai  Stegr«I“  >s  TTr  a"”  *”  “  lnpUt  “  the  integration.  The 
»•<««  iw  8  r,  ‘B»  ln  a  toPdown  sense,  above  the  derivative  commits- 
tion.  figure  6-2  shows  the  flow  of  derivative  computation. 


6. 1.2. 4  Output  Module 

This  module  produces  desired  output  describing  any  of  the  other  GO  Sub- 
system  functions.  The  output  may  be  printout,  plo£s,  or  d.M  s«f  me 

"X  j%the  “d“le  Is  the  0M  Is  dsta  set  and  the  primary 'oMput 

=  as r- 

research,  benchmark  and  contract  requirements  dictate. 


6.2  The  BUILD  Subsystem 


IJ\r1SnP^!K>aa  of  tha  BUIlD  Subsystem  is  to  prepare  an  executable  version 
SubsystM.  5'8'6"'  gUta  t~3  8hC"'"  the  functional  flow  of  the  BUILD 
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Figure  6-2.  Derivative  Compute  dons 
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6.2.1  BUILD  Data  Sets 

“li?.lnedyL“L^.‘eVJd  “pdateS  varlt"“  pr°8r“  libraries.  These 
maintained  as  partitioned  data  sets  on  direct  access  devices  It  flion 

These  "da  ta'sets'are ‘sequential”1  "T"  lihT&Tie8  and  CPP  da‘a  a«t. 

data  i!  !  a Ti  q  1  generally  a  tape.  Each  of  the  catalogued 

prefix  iual\f  ?  b  °W  ^C6Pt  thfe  Sy8te,n  Procedure  library  will  have  the 
are  IamLt  y  !?inaine  PGATWA6A,  e.g.,  FGATWA6A.  SOURCE.  All  data  sets 
are  ISrmanent.  All  except  AUCHIVE  are  pa  ^itioned  data  sets. 

fl‘  SS  \???tfln8  8°UrCe  COde  °f  a11  Pr°8ra“8  except 
OS/360  utilities  used  in  the  baseline  system. 

b.  OBJECT  -  Contains  object  modules  of  source  code  compiled 
from  SOURCE  members.  H 


c. 

d. 


h. 


LOAD  -  Contains  load  modules  of  OBJECT  members 

LEL  -  Contains  linkage  editor  code  for  forming  load 

modules  from  OBJECT  modules  and,  possibly,  other  load 
modules • 

JCL  -  Contains  Job  C  ntrol  Language  (JCL)  used  to 
execute  the  three  baseline  subsystems.  At  a  minimum, 
copies  of  the  procedures  for  executing  the  BUILD.  GO 
and  RESTORE  Subsystems  will  oe  members  of  this  data  set. 

ARCHIVE 1  -  Contains  a  copy  of  the  project  libraries 
resulting  from  a  pass  through  the  BUILD  Subsystem,  and 
a  copy  of  the  CPP  data  set  before  the  master  records 
are  updated  to  reflect  the  just  completed  BUILD  step. 

ARCHIVE2  -  Contains  a  copy  of  the  CPP  data  set  after 
the  master  records  have  been  updated  to  reflect  the 
BUILD  step. 

PROCLIB  -  The  operating  system  procedure  library.  JCL 
procedures  used  to  execute  jobs  must  be  copied  from 
project  libraries,  such  as  FGATWA6A. JCL,  to  PROCLIB 
before  they  can  be  executed. 


6.2.2  BUILD  Subsystem  Components 

8Jep  18  t0  read  and  Process  input  parameters.  These  parameters 
R°  rr  T®!’  executi°n  and  Precompile  parameters.  Execution  parame- 
"lriU8e  When  executing  the  GO  Subsystem.  Integration  time  step  is 
an  example  of  an  execution  parameter.  K 
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Precompile  parameters  can  be  subdivied  into  JCL  and  SOURCE  code  parameters . 
They  are  used  to  build  the  JCL  procedure  and  the  source  programs  that  will 
comprise  the  GO  Subsystem.  They  may  be  thought  of  as  execution  parameters 
for  the  BUILD  Subsystem.  JCL  parameters  are  used  to  compute  such  things 
as  data  set  record  length  and  blocksize.  Source  code  parameters  concern 
work  storage  sizes  end  segments  of  code  to  be  included  in  program  compila¬ 
tion.  Precompile  parameters  may  be  functions  of  other  precompile 
parameters . 

The  CPP  program  uses  input  parameters  to  update  the  CPP  master  and  data 
records,  and  checks  for  a  valid  and  complete  set  of  input.  If  the  data 
is  in  error,  CPP  sets  a  job  condition  code  of  99,  and  the  subsequent  job 
steps  are  not  executed.  The  CPP  program  also  prepares  a  temporary  data 
set  containing  utility  control  statements  to  be  used  for  archiving  the 
GO  Subsystem  once  it  is  built. 

At  the  completion  of  this  step,  CPP  data  records  containing  precompile 
parameters  have  been  created,  and  the  CPP  master  record  pointers  have 
been  modified  to  index  appropriate  data  records. 

The  next  step  is  to  assign  values  to  all  precompile  parameters  used  to 
build  the  GO  Subsystem  JCL  procedure  and  source  programs.  This  is  done 
in  a  PL/I  program  that  accesses  appropriate  CPP  data  records,  computes 
values  for  parameters  that  are  functions  of  other  parameters,  and  builds 
two  data  sets: 

a.  JCLPARM  -  A  member  of  the  project  JCL  library.  It 
contains  statements  assigning  values  to  all  JCL  pre- 
compile  parameters. 

b.  SOURCEPARM  -  A  member  of  the  project  SOURCE  library.  It 
contains  assignments  for  all  source  code  precompile 
parameters  used  in  the  Input,  Table  Generator,  Predictor, 
and  Output  modules. 

After  JCLPARM  and  SOURCE  PARM  are  written,  the  JCL  procedure  for  executing 
the  GO  Subsystem  must  be  built.  This  is  done  using  the  PL/I  precompiler 
facility.  A  program  executed  by  the  precompiler  has  the  following  features: 

a.  ZINCLUDE  JCL  Parameter  declarations. 

This  builds  in  line  code  identifying  all  JCL  parameters 
as  certain  precompiler  variable  typos 

b.  ZINCLUDE  JCLPARM. 

This  makes  the  contents  of  the  JCLPARM  data  set  built 
in  the  previous  step  inline  code.  Values  are  thus 
assigned  to  all  JCL  parameters  declared  via  the  pre¬ 
vious  ZINCLUDE. 
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Code  using  the  JCL  parameters  to  build  JCL 
statements. 


The  output  of  this  program  is  card  images  of  the  JCL  procedure  to  be  used 
for  executing  the  GO  Subsystem.  The  output  is  written  on  the  project  JCL 
library  as  member  WA6G0.  This  process  produces  a  JCL  procedure  with  data 
set  definitions  consistant  with  the  execution  and  precompile  parameters 
input  to  the  BUILD  Subsystem. 

Each  of  the  GO  Subsystem  program  modules  are  prepared  in  a  similar 
feshior. .  The  precompiler  is  used  for  a  first  pass  through  the  source  code 
to  resolve  references  to  precompiler  variables •  A  source  code  control 
segment  executed  by  the  precompiler  contains  ZINCLUDE  8  to. 

a.  Declare  source  code  parameters 


b.  Put  SOURCEPARM  contents  in  line 

c.  Define  work,  storage  areas  as  functions  of  precompile 
parameters 

d.  Define  I/O  data  sets  as  functions  of  precompile 
parame ters 


e.  Include  code  segments. 


The  code  segments  may  be  included  or  excluded  depending  on  values  of  pre¬ 
compile  parameters.  Further,  they  may  contain  expressions  using  precom¬ 
pile  variables.  The  output  of  the  precompiler  is  Source  code  resulting 
from  resolution  of  the  precompile  parameters  and  ready  for  compilation. 
This  code  is  compiled  and  object  modules  for  Input,  Table  Generator, 
Predictor  and  Output  are  written  on  the  project  OBJECT  librai”*/. 


The  object  modules  are  link  edited  and  executable  load  modules  written 
on  the  project  LOAD  library.  The  JCL  procedure  previously  built  is  copied 
to  the  operating  system  procedure  library  for  subsequent  execution. 

The  version  of  the  GO  subsystem  just  constructed  is  archived  for  possible 
future  restoration.  This  involves  copying  the  following  xnformation  onto 
a  sequential  data  set: 


a.  The  project  SOURCE  library 

b.  The  project  OBJECT  library 

c.  The  project  LOAD  library 

d.  The  project  JCL  library 

e.  The  CPP  data  set. 
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After  archiving,  the  BUILD  is  complete.  The  CPP  program  marks  and  updates 
the  master  records  accordingly,  and  the  updated  CPP  records  are  also 
archived. 


6.3  The  RESTORE  Subsystem 

The  RESTORE  Subsystem  retrieves  the  executing  portion  of  an  archived  GO 
Subsystem  and  restores  it  to  appropriate  libraries  for  subsequent  execu 
tion.  The  functional  flow  is  shown  in  Figure  6-4. 


6.3.1  RESTORE  Data  Sets 

The  RESTORE  Subsystem  accesses  permanent  data  sets  used  in  the  other  sub¬ 
systems,  and  builds  four  temporary  data  sets  for  its  own  use. 

a.  ARCHIVE 1  -  Sequential  data  set  containing  an  archived 
GO  Subsystem.  RESTORE  copies  the  JCL  and  LOAD  data  sets 
back  to  temporary  partitioned  direct  access  data  sets. 

b.  ARCHIVE2  -  Sequential  data  set  containing  CPP  data  set 
master  records.  RESTORE  copies  these  back  to  the  CPP 
data  set. 

c.  TEMPSYSIN  -  A  sequential  data  set  built  by  the  CPP 
program.  It  contains  utility  control  statements  used 
in  copying  from  ARCHIVEl  and  ARCHIVE 2  to  TEMPJCL , 

TEMPLOAD  and  TEMPCPP. 

d.  TEMPJCL  -  A  temporary  partitioned  data  set  containing 
the  JCL  library  copied  from  ARCHIVEl. 

e.  TEMPLOAD  -  A  temporary  partitioned  data  set  containing 
the  LOAD  library  copied  from  ARCHIVEl. 

f.  TEMPCPP  -  A  tf.nporary  direct  data  set  containing  the 
CPP  data  copied  from  ARCHIVE2. 

g.  PROCLIB  -  A  partitioned  data  set  containing  the  operating 
system  procedure  library.  The  GO  JCL  procedure  is  copied 
from  TEMPJCL  to  PROCLIB  so  that  it  can  be  executed. 

h.  LOAD  -  A  partitioned  data  set  containing  the  project  LOAD 
library.  The  Input,  Table  Generator,  Predictor,  and  Out¬ 
put  load  modules  are  copied  from  TEMPLOAD  to  LOAD  for 
subsequent  execution. 
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6.3.2  RESTORE  Subsystem  Components 

The  CPP  program  marks  the  CPP  system  master  record  "no-go"  to  indicate 
that  the  restore  process  is  not  complete  and  prevent  system  execution 
until  it  is.  CPP  also  builds  an  input  data  set  of  utility  control  cards 
for  retrieval  of  the  GO  Subsystem  in  the  next  step. 

The  next  step  recreates  the  LOAD  and  JCL  project  libraries  from  ARCHIVZl 
anu  the  CPP  data  set  from  ARCH1VE2. 

The  Input,  Table  Generation,  Predictor,  and  Output  load  modules  are 
copied  from  the  temporary,  restored  LOAD  data  set  to  the  project  LOAD 
library.  The  GO  JCL  procedure  is  copied  from  the  temporary  JCL  data  set 
to  the  system  procedure  library.  The  CPP  master  records  are  copied  from 
the  temporary  CPP  data  set  to  the  CPP  data  set. 

The  restored  GO  Subsystem  is  ready  for  execution.  The  CPP  program  marks 
the  system  master  record' accordingly,  and  indicates  that  archived  tables 
must  be  restored  when  the  GO  Subsystem  is  executed  for  the  first  time. 

This  table  restoration,  while  logically  part  of  the  RESTORE  step,  must  be 
done  in  the  GO  step.  The  JCL  (which  varies  from  version  to  version) 
describing  the  tables  is  contained  in  the  GO  procedure. 


6.4  Data  Set  Summary 

Table  6-1  enumerates  the  data  sets  used  in  the  baseline  system.  The  fol¬ 
lowing  classifications  are  used: 


a. 

Data 

Set  Organization 

1. 

Direct 

2. 

Partitioned 

3. 

Sequential 

b. 

Data 

Set  Type 

1. 

Permanent  -  exists  after  system 

executes 

2. 

Temporary  -  exists  only  during  subsystem  execution 

c. 

GO, 

BUILD,  RESTORE 

1. 

I  -  used  as  input  for  some  step 

in  the  subsystem 

2. 

0  -  output  produced  on  the  data 
in  the  subsystem. 

set  by  some  step 
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Table  6-] .  Data  Set  Summary 


Name 

Organization 

Type 

Go 

Build 

Restore 

CPP 

i" 

DIRECT 

PERM 

I/O 

I/O 

I/O 

«  B  fs 

O  ifl  U 
tl  b  it 

1  SOURCE 

PARTITIONED 

PERM 

I/O 

1-)  60  bJ 

0  0,0 

J  OBJECT 

PARTITIONED 

PERM 

o 

»-l  M  -W 

CM  CM  ij 

]  LOAD 

PARTITIONED 

PERM 

0 

LEL 

PARTITIONED 

PERM 

I 

_JCL 

PARTITIONED 

PERM 

I/O 

0) 

0) 

PR0CLIB 

PARTITIONED 

PERM 

o 

B  > 

pm, 

0)  *H 

■U  JS 
(0  u 

ARCHIVE1 

SEQUENTIAL 

PERM 

0 

I 

>.  n 
cn  <f 

ARCHIVE 2 

SEQUENTIAL 

PERM 

0 

I 

ARCHIVE3 

SEQUENTIAL 

PERM 

I/O 

m 

a 

RAW  DATA 

SEQUENTIAL (?) 

PERM 

I 

B  o 

flN  OUT 

0)  *H 

>  4J 

<0 
o 

W  iM 

SEQUENTIAL 

PERM 

I/O 

|  PM  IN 

SEQUENTIAL 

PERM 

I/O 

•33  s 

to  3  E 

trt  C 

PM  OUT 

SEQUENTIAL 

PERM 

I/O 

Q  o  o 

O  s  o 

ON  IN 

SEQUENTIAL 

PERM 

I/O 

ON  OUT 

SEQUENTIAL (?) 

PERM 

I/O 

TABLE  CONTROL 

SEQUENTIAL 

PERM 

0 

PREDINIT 

SEQUENTIAL 

PERM 

I/O 

u 

RADIAL  ONE 

DIRECT 

PERM 

I/O 

U 

4J 

U  O  CO 

RADIAL  TWO 

DIRECT 

PERM 

I/O 

•rl  H  CJ  J 

T3  3H  / 

0)  T3  ,£> 

PSLEGENDRE 

DIRECT 

PERM 

I/O 

MOW 

CM  2  H 

PSFOURIER 

DIRECT 

PERM 

I/O 

VSREAL 

DIRECT 

PERU 

I/O 

VSCOMPLEX 

DIRECT 

PERM 

I/O 

M  « 

O  U 

SPLEGENDRE 

DIRECT 

PERM 

I/O 

M  It  / 

O  t'J)  W  1 

•HO  I 

"3>!H  y 

P VALUE 

DIRECT 

TEMP 

I/O 

0>  K  ,£  \ 

bon 

PDIFF 

DIRECT 

TEMP 

I/O 

CM  S  H 

SVALUE 

DIRECT 

TEMP 

I/O 

c  DR 

DIRECT 

TEMP 

I/O 

<0 

u 

SDT 

DIRECT 

TEMP 

I/O 

<u  « 

M  Q 

O 

PDVAR 

DIRECT 

TEMP 

I/O 

Rest 

Work 

Sets 

s, 

TEMPJCL 

PARTITIONED 

TEMP 

I/O 

r 

TEMPLOAD 

PARTITIONED 

TEMP 

I/O 

fEMPCPP 

PARTITIONED 

TEMP 

I/O 

Hi 


Appendix  A.  PRECOMPILATION  CAPABILITIES  USING  A  PRECOMPILER 
IN  A  PROGRAM  PRODUCTION  ENVIRONMENT 

SE’S'.ssSs.*  srr1^8 a  pr°*r™ eha£  ^ »«« 

bilitv  DrovS  K  IhiS  18  similar  to  the  "macro  instruction"  capa¬ 
bilities  of  the  PL/rp^ecomenle(l)  ^  discussion  ls  United  to  capa- 
FORTRAN  to  allow  precompiling  in  thlt1£!U!£.bee"  lote8raee‘1  ”Uh 

The  capabilities  provided  by  this  package  allow: 

a-  Including  source  code  members  from  a  library 

b'  ofdff0'10"8  °f  arfay  di“ansions,  the  scope 

of  do-ioops  and  names  of  variables)  from  parameters 
set  at  compile  time. 

C.  A  "macro"  capability  whereby  the  user  can  extend  the 
linguage  syntactically  to  suit  hie  own  needs. 

d.  Conditional  compilation  of  code. 

With  these  facilities,  the  user  c»n 

a.  Eliminate  the  repetitious  task  of  duplicating  common 
code  (e.g.,  COMMON  statements,  DIMENSION  and  EQUIVALENT 

p  recoded  Members'!10'6  °°C‘’  by  S“ltabla  ■* 

b'  8Cjructurl"8  process  of  system  development  by 

including  dustily  code  sections  which  are  filled  in  as 
the  system  is  developed. 

C>  r"m'!”8  tlM  (whlle  Preserving  modularity)  by 

rh.J  dl?B  Cjde'  rather  than  calling  subroutines  with 
their  attendant  linkage  delay. 

d.  Compile  many  versions  of  a  program  by  setting  parame- 
ters.  without  altering  Source  code. 

®*  SfiJ?  "macJ08"  Horten  the  process  of  writing 
similar  code  sections.  ° 


a) 


"“T^Tug1!8^?  Intr°d“Ctl°”  “  tha  Compile-Time 
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In  the  following  discussion,  we  will  illustrate  these  capabilities  by 
means  of  examples.  In  order  to  avoid  an  irrelevent  discussion  of  the  pre¬ 
compiler’s  syntax,  the  preprocessor  statements  are  indicated  by  un  e 
lining  english  language  statements.  FORTRAN  will  be  used  as  the  source 
language,  although  the  facilities  are  also  part  of  PL/1  and  have  been 
extended  to  cover  OS/360  Job  Control  Language. 


Example  1.  Including  Code  from  a  source  library.  In  the  ,, 

example.  Source  library  LIBRARY  contains  a  member,  COHN,  with  the  System 
COMMON  statements.  The  programme?  avoids  having  to  duplicate  COMMON 
stater  ants  by  including  the  text  of  COMN: 


SOURCE  PROGRAM 
SUBROUTINE  SUB 
Include  COMN  from  LIBRARY 


MEMBER  COMN 
COMMON/A/X,  Y,  Z 
COMMON/B/U,  V,  W 


Additional  Code 


RETURN 

END 

The  text  produced  by  the  precompiler  is 
SUBROUTINE  SUB 
COMMON/A/X,  Y,  Z 
COMMON/B/U,  V,  W 

Additional  Code 
RETURN 


END 


Example  2.  Par ameti zing  array 
sion  of  array  A,  and  the  scope 
variable  "DIM" 


dimensions  and  other  variables.  The  dimen- 
of  the  do- loop,  are  from  a  preprocessor 
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SOURCE  PROGRAM 

!  |  !  1 

Set  Preprocessor  variable  DIM  ■  10 

,  '  i 

DIMENSION  A  (DIM) 

.  .  ‘  1  , 

DO  100  I  -  1,  DIM  1  i 

A  (I)  -  B  (I) 

100  CONTINUE 

The  text  produced  by  the  Precompiler  is 
DIMENSION  A  (10)  ,  i 

DO  100  I  -  1,  10 

i  1 

A  (I)  -  B  (I)  1 

100  CONTINUE  1  '  i  1 

Example  3.  Conditional  compilation  of  code.  If  the  preprocessor  variable 
DEBUG  is  set  to  PRODUCTION  the  code  to  call  PDUMP  is  to  be  skipped. 

i 

SOURCE  CODE  !  1  , 

Set  Preprocessor  variable  DEBUG  -  PRODUCTION  : 

SUBROUTINE  S(B) 

A-B  i  i 


If  Preprocessor  variable  DEBUG  **  PRODUCTION  then 

1 

1  1  I 

CALL  PDUMP  (A,B) 

Also 

>  i  .  1 

CALL  PRINT  (A) 

i 

END  ( 

Precompiler  Output 

1 

SUBROUTINE  S(B)  1 

A-B 

CALL  PRINT  (A) 

END  , 

A-3  '  ' 


I 


1 


| 


If  the  Preprocessor  variable  DEBUG  were  set  to,  e.g., 
"TEST"  in  the  first  line  above,  the  code  produced 
would  be: 

SUBROUTINE  S(B) 


A*'B 


CALL  PDUMP  (A,B) 

END 

Because  of  the  "else"  statement,  the  call  to  PRINT  was  skipped  and  the 
call  to  PDUMP  copied. 

Example  4.  Alternate  formulation  of  code.  In  this  example,  the  variable 
S  has  alternative  expressions.  Tn  one  case,  it  is  to  be  computed  as 
C*F*SIN (ALPHA),  in  another,  as  B*C/FUNC  (P).  This  is  accomplished  by 
setting  Preprocessor  variable  FORMULATION. 

SOURCE  CODE 


Set  Preprocessor  variable  FORMULATION  •  Fl 
If  Preprocessor  variable  FORMULATION  m  FI 
S  -  C*F*SIN  (ALPHA) 


Else 


S  -  B*C/FUNC  (P) 

OUTPUT  of  Precompiler  in  this  case  is 


S  -  C*F*SIN  (ALPHA) 

If  FORMULATION  were  set  to  some  other  value  the  other 
expression  would  be  copied  into  Source  code. 

Example  5.  A  macro  called  "STEP"  is  to  return  code  to  generate  the  com¬ 
monly  used  expression 

variable  ■■  *  variable  +  1. 
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SOURCE  CODE 


Macro  Definition  STEP(X):  Return  (X  ■  X  +  1) 

(This  line  states  that  the  macro  returns  the  form 
(argument  followed  by  sign  followed  by  argu¬ 
ment  followed  by  "+1"). 


STEP  (ALPHA) 
STEP  (BETA) 
STEP  (Z) 


END 

The  output  of  the  Precompiler  is 
ALPHA  -  ALPHA  +  1 
BETA  =  BETA  +  1 
Z  -  Z  +  1 


END 

The  programmer  is  spared  the  necessity  of  coding  each  statement  shown 
above.  In  this  sense  FORTRAN  has  been  extended  to  include  a  form  "STEP1 
among  its  keywords. 
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